Merge pull request #68 from kschan0214/dev1.2.2.5

Dev1.2.2.5

README.md

@@ -8,7 +8,7 @@
 
 The current GUI version is built to access the following toolboxes:
 - [MEDI (updated Jan 15, 2020)](http://weill.cornell.edu/mri/pages/qsm.html), 
-- [STI Suite (v3.0)](https://people.eecs.berkeley.edu/~chunlei.liu/software.html),
+- [STI Suite (v3.0)](https://chunleiliulab.github.io/software.html),
 - [FANSI (v3.0, released on 2021.10.15, i.e., commit b6ac1c9e)](https://gitlab.com/cmilovic/FANSI-toolbox/-/tree/b6ac1c9ea03380722ebe25a6dbef33fff4ea3700),  
 - [SEGUE](https://xip.uclb.com/i/software/SEGUE.html), and 
 - [nonlinear dipole inversion (NDI)](https://github.com/polakd/NDI_Toolbox),
@@ -41,27 +41,32 @@ If you have a more general question regarding the usage of SEPIA and/or other QS
 
 For full update log, please visit https://sepia-documentation.readthedocs.io/en/latest/getting_started/Release-note.html.
 
-### 1.2.2.4 (current master)
+### 1.2.2.5 (current master)
+* Fix the mismatch between SEPIA defined B0 direction and LPCNN when it is not along the z-direction
+* Fix the shared library issue when using ROMEO with latest versions of Matlab on Linux (see [here](https://github.com/korbinian90/ROMEO))
+* Allow user to define atlases' directory paths
+
+### 1.2.2.4 (commit 9083249)
 * Fix bug when importing SEPIA pipeline configuration files (sepia_config.m) to the GUI for using VSHARP and FANSI
 
-### 1.2.2.3 (current efde35b)
+### 1.2.2.3 (commit efde35b)
 * Fix bug when using BIDS compatible directory input where magnitude images did not utilise the rescale slope and intercept to obtain the true values for R2* mapping
 
-### 1.2.2.2 (current e53fd99)
+### 1.2.2.2 (commit e53fd99)
 * Fix bug when using BIDS compatible directory input where magnitude images did not utilise the rescale slope and intercept to obtain the true values for QSM
 
-### 1.2.2.1 (current 1f04298)
+### 1.2.2.1 (commit 1f04298)
 * Fix bug when using optimum weight total field computation with odd matrix size data
 
-### 1.2.2 (current d6bb60e)
+### 1.2.2 (commit d6bb60e)
 * Fix bug for non-double type input for MATLAB's strel function
 * Make sure all holes inside the ROI mask are filled after the background field removal step
 * ROI (brain) mask is applied on the fieldmap regardless of what method is chosen
 
-### 1.2.1.1 (current 941cd5b)
+### 1.2.1.1 (commit 941cd5b)
 * Enable option of GPU processing for FANSI and NDI
 
-### 1.2.1 (current 190dd44)
+### 1.2.1 (commit 190dd44)
 * Fix bug for data with odd-number matrix size
 * Fix bug for missing file when using R2* mapping with NLLS algorithm
 

SpecifyAtlasDirectory.m

@@ -0,0 +1,24 @@
+%% SpecifyAtlasDirectory.m
+%
+% Description: Specify directories of the atlas
+%
+% Kwok-Shing Chan @ MGH
+% [email protected]
+% Date created: 9 October 2023
+% Date modified: 
+%
+%% Specify the directories of the atlas here
+% please do not change
+SEPIA_ANALYSIS_SEGMENTATION_dir = fullfile(SEPIA_HOME,'analysis','segmentation');
+SEPIA_ATLAS_dir                 = fullfile(SEPIA_HOME,'atlas');
+
+% These are the default directories if you use the 'download_atlas.sh' shell script
+% If you prefer the atlases to be stored at a different paths, you may
+% modified the paths here
+
+% AHEAD atlas
+AHEAD_ATLAS_HOME                = fullfile(SEPIA_ATLAS_dir,'AHEAD_atlas');
+% CIT168 atlas
+CIT168_reinf_learn_ATLAS_HOME   = fullfile(SEPIA_ATLAS_dir,'CIT168_Reinf_Learn_v1.1.0');
+% MUSUS-100 atlas
+MuSus100_ATLAS_HOME             = fullfile(SEPIA_ATLAS_dir,'MuSus-100_Atlas');

addons/bfr/BFRnet/setup_BFRnet_environment.m

@@ -1,8 +1,8 @@
 %% This file specifies the Checkpoint of xQSM
 
 % Specify the directory of the xQSM code from Github
-deepMRI_HOME   	= '/project/3015069.05/bids/code/deepMRI/';
+deepMRI_HOME   	= '/home/common/matlab/sepia/external/deepMRI/';
 % If you have your own trained network, specify the file that contains the trained parameters
 % otherwise, uses the one provided with the tool
-checkpoints     = '/project/3015069.05/bids/code/deepMRI/BFRNet_data/checkpoints/BFRnet_L2_64PS_24BS_45Epo_NewHCmix.mat';
+checkpoints     = '/home/common/matlab/sepia/external/deepMRI/BFRNet_data/checkpoints/BFRnet_L2_64PS_24BS_45Epo_NewHCmix.mat';
 

addons/echo_combine/romeo/Wrapper_TotalField_ROMEO.m

@@ -40,6 +40,10 @@
 % parameters.output_dir = fullfile(tempdir, 'romeo_tmp');
 mkdir(parameters.output_dir);
 
+%% check shared library
+% 20231009: KC libraries with potential known issue
+if isunix; paths = getenv('LD_LIBRARY_PATH'); setenv('LD_LIBRARY_PATH'); end
+
 %% main
 [fieldmapUnwrapAllEchoes, totalField] = ROMEO(wrappedField, parameters);
 
@@ -70,6 +74,9 @@
 totalField                      = totalField *2*pi;
 totalField(isnan(totalField))   = 0;
 totalField(isinf(totalField))   = 0;
+
+% restore LD_LIBRARY_PATH environment
+if isunix; setenv('LD_LIBRARY_PATH', paths); end
 
 end
 

addons/phase_unwrap/romeo/Wrapper_Unwrap_ROMEO_singleEcho.m

@@ -42,12 +42,19 @@
 end
 mkdir(parameters.output_dir);
 
+%% check shared library
+% 20231009: KC libraries with potential known issue
+if isunix; paths = getenv('LD_LIBRARY_PATH'); setenv('LD_LIBRARY_PATH'); end
+
 %% main
 [unwrappedField, ~] = ROMEO(wrappedField, parameters);
 
 %% Delete tmp directory
 % Remove all temp output files and the temp folder
 rmdir(parameters.output_dir, 's')
+
+% restore LD_LIBRARY_PATH environment
+if isunix; setenv('LD_LIBRARY_PATH', paths); end
 
 end
 

addons/qsm/LPCNN/Wrapper_QSM_LPCNN.m

@@ -80,6 +80,8 @@
     mask_pad        = mask;
 end
 
+% 20231009: fix mismatch of B0dir between SEPIA and LPCNN
+b0dir = [b0dir(2), -b0dir(1), b0dir(3)];
 
 [C,~] = DipoleKernel(size(localField_pad),voxelSize,b0dir);
 

addons/qsm/LPCNN/setup_LPCNN_environment.m

@@ -3,9 +3,9 @@
 % dir_net             = fullfile(qsmnet_dir,'Checkpoints/');
 
 % Specify LPCNN_HOME corresponding to the code from Github
-LPCNN_HOME          = '/project/3015069.05/bids/code/LPCNN/';
+LPCNN_HOME          = '/home/common/matlab/sepia/external/LPCNN/LPCNN_v1.0';
 % Specify the Python environment that has QSMnet+ installed
-python_interpreter  = '/project/3015069.05/bids/code/LPCNN/lpcnn-env/bin/python';
+python_interpreter  = '/home/common/matlab/sepia/external/LPCNN/LPCNN_v1.0/lpcnn-env/bin/python';
 % Specify the directory that contains the training parameters
 checkpoint_fn       = fullfile(LPCNN_HOME,'checkpoints','lpcnn_test_Bmodel.pkl');
 % checkpoint_fn       = '/project/3015069.05/bids/code/LPCNN/checkpoints/lpcnn_test_Bmodel.pkl';

addons/qsm/QSMnet/setup_qsmnet_environment.m

@@ -1,8 +1,8 @@
 %% This file specifies the Python environment and Checkpoint of QSMnet+
 
 % Specify the directory of the QSMnet+ code from Github
-QSMnet_HOME         = '/project/3015069.05/bids/code/QSMnet/';
+QSMnet_HOME         = '/home/common/matlab/sepia/external/QSMnet/QSMnet-9ca283a';
 % Specify the Python environment that has QSMnet+ installed
-python_interpreter  = '/project/3015069.05/bids/code/QSMnet/qsmnet/bin/python';
+python_interpreter  = '/home/common/matlab/sepia/external/QSMnet/QSMnet-9ca283a/qsmnet/bin/python';
 % Specify the directory that contains the trained parameters
 dir_net             = fullfile(QSMnet_HOME,'Checkpoints/');

addons/qsm/xQSM/setup_xQSM_environment.m

@@ -1,7 +1,7 @@
 %% This file specifies the Checkpoint of xQSM
 
 % Specify the directory of the xQSM code from Github
-deepMRI_HOME   	= '/project/3015069.05/bids/code/deepMRI/';
+deepMRI_HOME   	= '/home/common/matlab/sepia/external/deepMRI/';
 % If you have your own trained network, specify the file that contains the trained parameters
 % otherwise, uses the one provided with the tool
 checkpoint_dir  = fullfile(deepMRI_HOME, 'xQSM_data','checkpoints');

analysis/segmentation/get_AHEAD_atlas_labels.m

@@ -18,7 +18,7 @@
 % Kwok-shing Chan @ DCCN
 % [email protected]
 % Date created: 6 October 2022
-% Date last modified:
+% Date modified: 9 October 2023
 %
 %
 function get_AHEAD_atlas_labels(input,output_dir,algorParam)
@@ -36,9 +36,11 @@ function get_AHEAD_atlas_labels(input,output_dir,algorParam)
 if status == 127
     setenv('PATH', [getenv('PATH') ':' ANTS_HOME]);
 end
-SEPIA_ANALYSIS_SEGMENTATION_dir = fullfile(SEPIA_HOME,'analysis','segmentation');
-SEPIA_ATLAS_dir                 = fullfile(SEPIA_HOME,'atlas');
-AHEAD_ATLAS_dir                 = fullfile(SEPIA_ATLAS_dir,'AHEAD_atlas');
+% get atlas directory
+SpecifyAtlasDirectory;
+% SEPIA_ANALYSIS_SEGMENTATION_dir = fullfile(SEPIA_HOME,'analysis','segmentation');
+% SEPIA_ATLAS_dir                 = fullfile(SEPIA_HOME,'atlas');
+% AHEAD_ATLAS_HOME                = fullfile(SEPIA_ATLAS_dir,'AHEAD_atlas');
 
 output_tmp_dir = fullfile(output_dir,'AHEAD_intermediate_files',filesep);
 if ~exist(output_dir,'dir')
@@ -81,11 +83,11 @@ function get_AHEAD_atlas_labels(input,output_dir,algorParam)
         save_nii_img_only(T1w_nii,img_hybrid_nii, img_hybrid );
 
         % create hybrid T1w+QSM image on the atlas data
-        t1w         = load_nii_img_only(fullfile(AHEAD_ATLAS_dir,'Templates_mni09b','ahead_final_med_r1map_n104_mni09b.nii.gz'));
-        Chi_t1w     = load_nii_img_only(fullfile(AHEAD_ATLAS_dir,'Templates_mni09b','ahead_final_med_qsm_n104_mni09b.nii.gz'));
+        t1w         = load_nii_img_only(fullfile(AHEAD_ATLAS_HOME,'Templates_mni09b','ahead_final_med_r1map_n104_mni09b.nii.gz'));
+        Chi_t1w     = load_nii_img_only(fullfile(AHEAD_ATLAS_HOME,'Templates_mni09b','ahead_final_med_qsm_n104_mni09b.nii.gz'));
         mask_t1     = t1w ~= 0;
         img_hybrid  = compute_hybrid_t1w_chi(t1w, Chi_t1w, mask_t1, 400, 0.5 );
-        save_nii_img_only(fullfile(AHEAD_ATLAS_dir,'Templates_mni09b','ahead_final_med_r1map_n104_mni09b.nii.gz'),template_hybrid_nii, img_hybrid);
+        save_nii_img_only(fullfile(AHEAD_ATLAS_HOME,'Templates_mni09b','ahead_final_med_r1map_n104_mni09b.nii.gz'),template_hybrid_nii, img_hybrid);
         clear img_hybrid Chi_t1w t1w mask_t1
 
         % Step 4: register hybrid image to template space, non linear transform
@@ -99,20 +101,20 @@ function get_AHEAD_atlas_labels(input,output_dir,algorParam)
         t1_2_t1wTemplate_mat                = fullfile(output_tmp_dir, 'T1w_2_Template_hybrid_image_0GenericAffine.mat');
         t1_2_t1wTemplate_inverseWrap_nii    = fullfile(output_tmp_dir, 'T1w_2_Template_hybrid_image_1InverseWarp.nii.gz');
 
-        mask_list = dir(fullfile(AHEAD_ATLAS_dir,'structures_mni09b','*_mask-*nii*'));
+        mask_list = dir(fullfile(AHEAD_ATLAS_HOME,'structures_mni09b','*_mask-*nii*'));
         mode_interp = 1; % 1=GenericLabel; 2=linear
         mode_4D     = 0; % 0=3D; 3=4D; equivalent to option -e in antsApplyTransforms see ANTs doc
         for k = 1:numel(mask_list)
-            label_nii   = fullfile(AHEAD_ATLAS_dir,'structures_mni09b',mask_list(k).name);
+            label_nii   = fullfile(AHEAD_ATLAS_HOME,'structures_mni09b',mask_list(k).name);
             cmd = ['sh ' shell_script ' ' output_tmp_dir ' ' num2str(mode_interp) ' ' num2str(mode_4D) ' ' label_nii ' ' Chi_nii ' ' gre_2_T1w_mat ' ' t1_2_t1wTemplate_mat ' ' t1_2_t1wTemplate_inverseWrap_nii];
             system(cmd);
         end
 
-        prob_list = dir(fullfile(AHEAD_ATLAS_dir,'structures_mni09b','*_proba-*nii*'));
+        prob_list = dir(fullfile(AHEAD_ATLAS_HOME,'structures_mni09b','*_proba-*nii*'));
         mode_interp = 2; % 1=GenericLabel; 2=linear
         mode_4D     = 0; % 0=3D; 3=4D; equivalent to option -e in antsApplyTransforms see ANTs doc
         for k = 1:numel(prob_list)
-            label_nii   = fullfile(AHEAD_ATLAS_dir,'structures_mni09b',prob_list(k).name);
+            label_nii   = fullfile(AHEAD_ATLAS_HOME,'structures_mni09b',prob_list(k).name);
             cmd = ['sh ' shell_script ' ' output_tmp_dir ' ' num2str(mode_interp) ' ' num2str(mode_4D) ' ' label_nii ' ' Chi_nii ' ' gre_2_T1w_mat ' ' t1_2_t1wTemplate_mat ' ' t1_2_t1wTemplate_inverseWrap_nii];
             system(cmd);
         end
@@ -126,20 +128,20 @@ function get_AHEAD_atlas_labels(input,output_dir,algorParam)
         % Step 1: Apply tranformation
         shell_script = fullfile(SEPIA_ANALYSIS_SEGMENTATION_dir,'ANTs_gre_2_t1wAtlas_applyTransform.sh');
 
-        mask_list = dir(fullfile(AHEAD_ATLAS_dir,'structures_mni09b','*_mask-*nii*'));
+        mask_list = dir(fullfile(AHEAD_ATLAS_HOME,'structures_mni09b','*_mask-*nii*'));
         mode_interp = 1; % 1=GenericLabel; 2=linear
         mode_4D     = 0; % 0=3D; 3=4D; equivalent to option -e in antsApplyTransforms see ANTs doc
         for k = 1:numel(mask_list)
-            label_nii   = fullfile(AHEAD_ATLAS_dir,'structures_mni09b',mask_list(k).name);
+            label_nii   = fullfile(AHEAD_ATLAS_HOME,'structures_mni09b',mask_list(k).name);
             cmd = ['sh ' shell_script ' ' output_tmp_dir ' ' num2str(mode_interp) ' ' num2str(mode_4D) ' ' label_nii ' ' GRE_nii ' ' gre_2_T1w_mat ' ' t1_2_t1wTemplate_mat ' ' t1_2_t1wTemplate_inverseWrap_nii];
             system(cmd);
         end
 
-        prob_list = dir(fullfile(AHEAD_ATLAS_dir,'structures_mni09b','*_proba-*nii*'));
+        prob_list = dir(fullfile(AHEAD_ATLAS_HOME,'structures_mni09b','*_proba-*nii*'));
         mode_4D     = 0; % 0=3D; 3=4D; equivalent to option -e in antsApplyTransforms see ANTs doc
         mode_interp = 2; % 1=GenericLabel; 2=linear
         for k = 1:numel(prob_list)
-            label_nii   = fullfile(AHEAD_ATLAS_dir,'structures_mni09b',prob_list(k).name);
+            label_nii   = fullfile(AHEAD_ATLAS_HOME,'structures_mni09b',prob_list(k).name);
             cmd = ['sh ' shell_script ' ' output_tmp_dir ' ' num2str(mode_interp) ' ' num2str(mode_4D) ' ' label_nii ' ' GRE_nii ' ' gre_2_T1w_mat ' ' t1_2_t1wTemplate_mat ' ' t1_2_t1wTemplate_inverseWrap_nii];
             system(cmd);
         end

analysis/segmentation/get_CIT168_reinf_learn_labels.m

@@ -18,7 +18,7 @@
 % Kwok-shing Chan @ DCCN
 % [email protected]
 % Date created: 6 October 2022
-% Date last modified:
+% Date modified: 9 October 2023
 %
 %
 function get_CIT168_reinf_learn_labels(input,output_dir,algorParam)
@@ -36,11 +36,13 @@ function get_CIT168_reinf_learn_labels(input,output_dir,algorParam)
 if status == 127
     setenv('PATH', [getenv('PATH') ':' ANTS_HOME]);
 end
-SEPIA_ANALYSIS_SEGMENTATION_dir = fullfile(SEPIA_HOME,'analysis','segmentation');
-SEPIA_ATLAS_dir                 = fullfile(SEPIA_HOME,'atlas');
-CIT168_reinf_learn_ATLAS_dir    = fullfile(SEPIA_ATLAS_dir,'CIT168_Reinf_Learn_v1.1.0');
+% get atlas directory
+SpecifyAtlasDirectory;
+% SEPIA_ANALYSIS_SEGMENTATION_dir = fullfile(SEPIA_HOME,'analysis','segmentation');
+% SEPIA_ATLAS_dir                 = fullfile(SEPIA_HOME,'atlas');
+% CIT168_reinf_learn_ATLAS_HOME   = fullfile(SEPIA_ATLAS_dir,'CIT168_Reinf_Learn_v1.1.0');
 
-template_nii = fullfile(CIT168_reinf_learn_ATLAS_dir,'MNI152-Nonlin-Asym-2009c','CIT168toMNI152-2009c_T1w_brain.nii.gz');
+template_nii = fullfile(CIT168_reinf_learn_ATLAS_HOME,'MNI152-Nonlin-Asym-2009c','CIT168toMNI152-2009c_T1w_brain.nii.gz');
 
 output_tmp_dir = fullfile(output_dir,'CIT168_reinf_learn_intermediate_files',filesep);
 if ~exist(output_dir,'dir')
@@ -75,14 +77,14 @@ function get_CIT168_reinf_learn_labels(input,output_dir,algorParam)
         % deterministic labels
         mode_interp = 1; % 1=GenericLabel; 2=linear
         mode_4D     = 0; % 0=3D; 3=4D; equivalent to option -e in antsApplyTransforms see ANTs doc
-        label_nii   = fullfile(CIT168_reinf_learn_ATLAS_dir,'MNI152-Nonlin-Asym-2009c','CIT168toMNI152-2009c_det.nii.gz');
+        label_nii   = fullfile(CIT168_reinf_learn_ATLAS_HOME,'MNI152-Nonlin-Asym-2009c','CIT168toMNI152-2009c_det.nii.gz');
         cmd = ['sh ' shell_script ' ' output_tmp_dir ' ' num2str(mode_interp) ' ' num2str(mode_4D) ' ' label_nii ' ' GRE_nii ' ' gre_2_T1w_mat ' ' t1_2_t1wTemplate_mat ' ' t1_2_t1wTemplate_inverseWrap_nii];
         system(cmd);
 
         % probabilistic labels
         mode_interp = 2; % 1=GenericLabel; 2=linear
         mode_4D     = 3; % 0=3D; 3=4D; equivalent to option -e in antsApplyTransforms see ANTs doc
-        label_nii   = fullfile(CIT168_reinf_learn_ATLAS_dir,'MNI152-Nonlin-Asym-2009c','CIT168toMNI152-2009c_prob.nii.gz');
+        label_nii   = fullfile(CIT168_reinf_learn_ATLAS_HOME,'MNI152-Nonlin-Asym-2009c','CIT168toMNI152-2009c_prob.nii.gz');
         cmd = ['sh ' shell_script ' ' output_tmp_dir ' ' num2str(mode_interp) ' ' num2str(mode_4D) ' ' label_nii ' ' GRE_nii ' ' gre_2_T1w_mat ' ' t1_2_t1wTemplate_mat ' ' t1_2_t1wTemplate_inverseWrap_nii];
         system(cmd);
 
@@ -97,22 +99,22 @@ function get_CIT168_reinf_learn_labels(input,output_dir,algorParam)
         % deterministic labels
         mode_interp = 1; % 1=GenericLabel; 2=linear
         mode_4D     = 0; % 0=3D; 3=4D; equivalent to option -e in antsApplyTransforms see ANTs doc
-        label_nii   = fullfile(CIT168_reinf_learn_ATLAS_dir,'MNI152-Nonlin-Asym-2009c','CIT168toMNI152-2009c_det.nii.gz');
+        label_nii   = fullfile(CIT168_reinf_learn_ATLAS_HOME,'MNI152-Nonlin-Asym-2009c','CIT168toMNI152-2009c_det.nii.gz');
         cmd = ['sh ' shell_script ' ' output_tmp_dir ' ' num2str(mode_interp) ' ' num2str(mode_4D) ' ' label_nii ' ' GRE_nii ' ' gre_2_T1w_mat ' ' t1_2_t1wTemplate_mat ' ' t1_2_t1wTemplate_inverseWrap_nii];
         system(cmd);
 
         % probabilistic labels
         mode_interp = 2; % 1=GenericLabel; 2=linear
         mode_4D     = 3; % 0=3D; 3=4D; equivalent to option -e in antsApplyTransforms see ANTs doc
-        label_nii   = fullfile(CIT168_reinf_learn_ATLAS_dir,'MNI152-Nonlin-Asym-2009c','CIT168toMNI152-2009c_prob.nii.gz');
+        label_nii   = fullfile(CIT168_reinf_learn_ATLAS_HOME,'MNI152-Nonlin-Asym-2009c','CIT168toMNI152-2009c_prob.nii.gz');
         cmd = ['sh ' shell_script ' ' output_tmp_dir ' ' num2str(mode_interp) ' ' num2str(mode_4D) ' ' label_nii ' ' GRE_nii ' ' gre_2_T1w_mat ' ' t1_2_t1wTemplate_mat ' ' t1_2_t1wTemplate_inverseWrap_nii];
         system(cmd);
 
 end
 
 copyfile(fullfile(output_tmp_dir,'CIT168toMNI152-2009c_det_2gre.nii.gz'),fullfile(output_dir));
 copyfile(fullfile(output_tmp_dir,'CIT168toMNI152-2009c_prob_2gre.nii.gz'),fullfile(output_dir));
-copyfile(fullfile(CIT168_reinf_learn_ATLAS_dir,'labels.txt'),fullfile(output_dir,'CIT168_reinf_learn_labels.txt'));
+copyfile(fullfile(CIT168_reinf_learn_ATLAS_HOME,'labels.txt'),fullfile(output_dir,'CIT168_reinf_learn_labels.txt'));
 
 disp('Segmentation is done!')