diff --git a/Changelog.md b/Changelog.md new file mode 100644 index 0000000..a760073 --- /dev/null +++ b/Changelog.md @@ -0,0 +1,59 @@ +New in this version [1.7]: +- Initial support for calculating structural connectivity (based on tools from EDTI) + +Update 11-12-2024: +- Bug fixes +- Updated compatibility to the latest version of WhiteMatterAnalysis +- Initial support for tractography with Scilpy tools + +Update 17-02-2023: +- This is the last stable release before a major change in structure. From the next release, all main functionalities will be come available as standalone commands and be designed to fully rely on Niftis only. +- Fixes to the handling of Nifti headers +- Automated brain extraction using image registration +- added new class MRT_Library that will collect all own methods currently in MRTTrack and MRTQuant (transition in progress). This will clean MRTTrack and MRTQuant which are meant as documented interfaces. +- Dramatic speed up of GRL and mFOD (major revision of MRTTrack.PerformDeconv) +- Changes to MRTQuant.LoadNifti and SaveNifti, which now allows to preserve the original header of NIFTIs (if required). In future releases, no modifications of the headers will be required. +- Update of Neuro.m to support the latest version of Elastix +- Compressed niftis are now (un)compressed in a temporary folder +- Added functions to track multiple FODs simultaneously +- Changes to the automatic handling of Q-S form of NIFTIs +- New mFOD fiber tracking +- New TerminateTractsWithFraction approach with 1 voxel tolerance +- Moved Trackers to aanother location +- Fixes to SHPrecomp +- Fixes to MRIToolkitInit +- Tractography can now use parallel computing +- Handling of NIFTI scaling factors +- Failsafe mode for registration +- Support for FSL Eddy +- Support for FOD normalization +- One call command to run GRL +- Fix for MSCSD (probably still needs amendment on the detection of data shells) +- Fix fot CAT12 processing +- command line fixes +- standardized preproc includes conformation +- fixes to conformation +- support for offset in DKI for GRL +- fixes to DKI export + + +Update 09-2021: +- [NEW] Checkout the [tutorial](https://github.com/delucaal/MRIToolkit/tree/master/GettingStarted/Diffusion/7_WMA) on how to install/use automated tractography clustering (WMA) algorithm with CSD/GRL/mFOD +- Fixes to MSCSD-related functions +- Fixes to WM automated clustering +- Support for Q-form in MRTTrack.ConformSpatialDimensions() +- New FOD scaling mode for GRL/mFOD +- Renaming of some command line functions +- Added an option to customise the number of dRL iterations +- New functions to estimate SNR from b=0s/mm2 images (MRTQuant) + + +Update 17-02-2021: +- Run the CAT12 automatic pipeline for T1 images +- Code re-organization into two main classes: MRTQuant (preprocessing/DTI/DKI) and MRTrack (Tractography related) +- Support for automatic fiber clustering (see below for reference) +- Early support for integration with Python (needed for the point above) +- Support for VTK poly data 4.2 and 5.1 (also needed for the clustering) +- Added a robust option to GRL/mFOD deconvolution +- Initial support for storing the NIFTI Q/S form (to improve interoperability with other tools, not implemented yet) +- Integration with CAT12 \ No newline at end of file diff --git a/Core/MRT_Library.m b/Core/MRT_Library.m index b1cfed2..06097e0 100644 --- a/Core/MRT_Library.m +++ b/Core/MRT_Library.m @@ -2065,6 +2065,20 @@ function WriteElastixParameters(fout,pars) end end + % Make an atlas file EDTI compatible + function AdaptAtlasEDTIStyle(parc_file,labels,ids,output) + pf = MRTQuant.LoadNifti(parc_file); + pf_new = pf; + pf_new.img = zeros(size(pf_new.img)); + labels_txt = fopen([output '_labels.txt'],'wt'); + for ij=1:length(ids) + pf_new.img(pf.img == ids(ij)) = ij; + fprintf(labels_txt,'1 %s%s',strrep(labels{ij},' ','_'),newline); + end + fclose(labels_txt); + MRTQuant.WriteNifti(pf_new,[output '_atlas.nii']); + end + % To be completed function fs_lut = FreeSurferLabelsInParcFile(parc_file) fs_lut = []; diff --git a/Core/mk_init.m b/Core/mk_init.m index 47388c2..aac1065 100644 --- a/Core/mk_init.m +++ b/Core/mk_init.m @@ -10,8 +10,9 @@ % 24/10/2020: creation - v1.1 % 01/01/2021: update v1.2 % 12/01/2024: update v1.3 +% 18/12/2024: update v1.4 global MRIToolkit; -MRIToolkit.core_version = 1.3; +MRIToolkit.core_version = 1.4; addpath(get_executed_file_path()) addpath(fullfile(get_executed_file_path(),'Trackers')) \ No newline at end of file diff --git a/ExploreDTIInterface/EDTI_Library.m b/ExploreDTIInterface/EDTI_Library.m index cb14088..471b89d 100644 --- a/ExploreDTIInterface/EDTI_Library.m +++ b/ExploreDTIInterface/EDTI_Library.m @@ -4615,6 +4615,8 @@ function E_DTI_model_fit(f_DWI,f_BM,fnam,Mask_par,perm,flip,fit_mode, dki_fit, d try hdr = load_untouch_header_only(f_DWI); catch + pf = dir([strrep(strrep(f_DWI,'.nii',''),'.gz','') '.nii*']); + hdr = load_untouch_header_only(fullfile(pf.folder,pf.name)); end b = textread(f_BM); @@ -8429,6 +8431,882 @@ function E_DTI_Convert_mat_2_nii(fin,folout,LS) seedPoint = seedPoint(:,1:i); end + % From ExploreDTI: network connectivity analysis based on tractography and + % brain parcellations in native space + function E_DTI_Network_analysis_from_ROI_L(files,tract_files,label_files,labels_descriptor,export_tracts,pass_or_end_or_both,out_folder,parallel_computing) + d = pwd; + try + [Lab_B, Lab] = textread(labels_descriptor,'%f%s','delimiter','\n'); + catch + disp(['Error reading label names file: ' full_L_n3]) + cd(d) + return; + end + + if ~all(any([Lab_B==1 Lab_B==0],2)) + error(['Error reading label names file: ' labels_descriptor]) + end + + if sum(Lab_B)<3 + error(['There should be at least 3 non-zero labels in file: ' labels_descriptor]) + end + + if (export_tracts == false) + selected_labels = []; + else + selected_labels = EDTI_Library.E_DTI_Get_NA_Tract_selection(Lab(Lab_B==1)); + end + + if(strcmpi(pass_or_end_or_both,'pass')) + ACh = 1; + elseif (strcmpi(pass_or_end_or_both,'end')) + ACh = 2; + elseif(strcmpi(pass_or_end_or_both,'both')) + ACh = 3; + end + + %ADL + % files = files(complete_files); + % label_files = label_files(complete_files); + vol_files = cell(size(files)); %this could be used to extract values of map from another modality for a specific connection %vol_files(complete_files); + vol_ext = 'nii'; + + if(parallel_computing==0) + + disp(' ') + tc = 0; + + for i=1:length(files) + + tic; + + disp('Processing file: ') + disp(files{i}) + + try + suc = EDTI_Library.E_DTI_Network_analysis_from_ROI_L_exe(files{i},tract_files{i},... + vol_files{i}, label_files{i}, out_folder, selected_labels, ACh, Lab, Lab_B, vol_ext); + catch me + disp('Unexpected error for file:') + disp(files{i}) + disp(me.message) + end + + if suc==0 + disp('Processing not successful.') + else + disp('Processing successful.') + end + + t=toc; + m=t/60; + disp(['Computation time was ' num2str(m) ' min.']) + + tc = tc + m; + + disp(' ') + + end + + disp('Network analysis completed.') + + disp(['Total computation time was ' num2str(tc) ' min.']) + + cd(d) + + + else + + disp(' ') + + tic; + + parfor i=1:length(files) + + disp('Processing file: ') + disp(files{i}) + + try + EDTI_Library.E_DTI_Network_analysis_from_ROI_L_exe(files{i},tract_files{i},... + vol_files{i}, label_files{i}, out_folder, selected_labels, ACh, Lab, Lab_B, vol_ext); + catch me + disp('Unexpected error for file:') + disp(files{i}) + disp(me.message) + end + + % if suc==0 + % disp('Processing not successful.') + % else + % disp('Processing successful.') + % end + + end + + t=toc; + m=t/60; + disp(['Computation time was ' num2str(m) ' min.']) + + disp('Network analysis completed.') + + cd(d) + + end + end + + % From ExploreDTI: helper function for connectivity analysis + function s = E_DTI_Get_NA_Tract_selection(L) + + [L_s,ind] = sortrows(L); + + cn=0; + + for i=1:length(L_s) + for j=1:length(L_s) + if i==j + nL = L_s{i}; + else + nL = [L_s{i} '_conn_' L_s{j}]; + end + + cn = cn+1; + Lab{cn} = nL; + + end + end + + s = 1:length(Lab); + [I, J] = ind2sub([length(L) length(L)],s); + + s = [J(:) I(:)]; + + for i=1:size(s,1) + s(i,:) = [ind(s(i,1)) ind(s(i,2))]; + end + end + + % From ExploreDTI: Helper function for connectivity analysis + function suc = E_DTI_Network_analysis_from_ROI_L_exe(dti_f,tract_f,... + vol_f, label_f, out_f, selected_labels, ACh, Lab, Lab_B, vol_ext) + + L = Lab(Lab_B==1); + + fiL = find(Lab_B); + + suc = 1; + + if isempty(selected_labels) + save_o = 0; + else + save_o = 1; + end + + if exist(tract_f,'file')~=2 + suc = 0; + disp(['File: ' tract_f ' does not exist.']) + return; + end + + try + [A_L_trafo, VDims] = EDTI_Library.E_DTI_read_nifti_file(label_f); + catch + suc = 0; + disp(['Problems reading file: ' label_f ' ...']) + return; + end + + for i=1:length(L) + + mask{i} = repmat(0>1,size(A_L_trafo)); + mask{i}(A_L_trafo==fiL(i))=1; + + end + + disp('Analyzing labels...') + if ACh==1 || ACh==2 + F_List = EDTI_Library.E_DTI_Get_F_List_from_tracts_and_masks(tract_f,mask,ACh); + else + F_List_1 = EDTI_Library.E_DTI_Get_F_List_from_tracts_and_masks(tract_f,mask,1); + disp(' ') + F_List_2 = EDTI_Library.E_DTI_Get_F_List_from_tracts_and_masks(tract_f,mask,2); + end + disp('Done analyzing labels.') + + disp(' ') + + if save_o==1 + disp('Saving fiber bundles...') + if ACh==1 || ACh==2 + EDTI_Library.E_DTI_save_Network_analysis_temp_tracts(tract_f,F_List,out_f,L,selected_labels,ACh); + else + EDTI_Library.E_DTI_save_Network_analysis_temp_tracts(tract_f,F_List_1,out_f,L,selected_labels,1); + disp(' ') + EDTI_Library.E_DTI_save_Network_analysis_temp_tracts(tract_f,F_List_2,out_f,L,selected_labels,2); + end + disp('Done saving fiber bundles.') + end + + disp(' ') + + disp('Calculating metrics...') + if ACh==1 || ACh==2 + [Conn_matrices, Conn_Lab] = ... + EDTI_Library.E_DTI_Get_Network_analysis_conn_mat_1(tract_f,F_List,L,vol_f,dti_f,vol_ext,ACh,mask); + else + [Conn_matrices_1, Conn_Lab_1] = ... + EDTI_Library.E_DTI_Get_Network_analysis_conn_mat_1(tract_f,F_List_1,L,vol_f,dti_f,vol_ext,1,mask); + disp(' ') + [Conn_matrices_2, Conn_Lab_2] = ... + EDTI_Library.E_DTI_Get_Network_analysis_conn_mat_1(tract_f,F_List_2,L,vol_f,dti_f,vol_ext,2,mask); + end + disp('Done calculating metrics.') + + disp(' ') + + if ACh==1 || ACh==2 + Out_names = EDTI_Library.E_DTI_get_NA_final_output_names(out_f,tract_f,Conn_Lab); + else + Out_names_1 = EDTI_Library.E_DTI_get_NA_final_output_names(out_f,tract_f,Conn_Lab_1); + Out_names_2 = EDTI_Library.E_DTI_get_NA_final_output_names(out_f,tract_f,Conn_Lab_2); + end + + if ACh==1 || ACh==2 + for i=1:length(Out_names) + CM = Conn_matrices{i}; + save(Out_names{i},'CM') + end + else + for i=1:length(Out_names_1) + CM = Conn_matrices_1{i}; + save(Out_names_1{i},'CM') + CM = Conn_matrices_2{i}; + save(Out_names_2{i},'CM') + end + end + + CM = repmat(nan,[length(mask) length(mask)]); + CM = EDTI_Library.E_DTI_Get_IL_dist(mask); + + for i=1:length(mask) + for j=i:length(mask) + CM(j,i) = CM(i,j); + end + end + + CM_ext{1} = '_inter_label_distances'; + Out_name = EDTI_Library.E_DTI_get_NA_final_output_names(out_f,tract_f,CM_ext); + save(Out_name{1},'CM') + end + + % From ExploreDTI: Helper function for connectivity analysis + function ndx = E_DTI_sub2ind(siz,a,b,c) + %Compute linear indices + k = [1 cumprod(siz(1:end-1))]; + ndx = 1; + + ndx = ndx + (a-1)*k(1); + ndx = ndx + (b-1)*k(2); + ndx = ndx + (c-1)*k(3); + + end + + % From ExploreDTI: Helper function for connectivity analysis + function Fin_List = E_DTI_Get_F_List_from_tracts_and_masks(tract_f,AND_mask,ACh) + + load(tract_f,'Tracts','VDims','TractMask'); + MDims = size(TractMask); + + Lt = length(Tracts); + Lts = zeros(Lt,1); + + if ACh==2 + for i = 1:Lt + Tracts{i} = Tracts{i}([1 end],:); + end + end + + for i = 1:Lt + Lts(i,1)=size(Tracts{i},1); + end + CLts = [0 ;cumsum(Lts)]; + + Tmatx = repmat(single(0),[sum(Lts) 1]); + Tmaty = repmat(single(0),[sum(Lts) 1]); + Tmatz = repmat(single(0),[sum(Lts) 1]); + + for i=1:Lt + Tmatx(CLts(i)+1:CLts(i+1),1) = Tracts{i}(:,1); + Tmaty(CLts(i)+1:CLts(i+1),1) = Tracts{i}(:,2); + Tmatz(CLts(i)+1:CLts(i+1),1) = Tracts{i}(:,3); + end + + Tmatx = Tmatx/VDims(1); + Tmatx = round(Tmatx); + Tmaty = Tmaty/VDims(2); + Tmaty = round(Tmaty); + Tmatz = Tmatz/VDims(3); + Tmatz = round(Tmatz); + + TInd = EDTI_Library.E_DTI_sub2ind(MDims,Tmatx,Tmaty,Tmatz); + clear Tmatx Tmaty Tmatz; + + + for i=1:length(AND_mask) + if ACh==1 + % disp(['Analyzing label ' num2str(i) ' of ' num2str(length(AND_mask)) ' (PASS)...']) + else + % disp(['Analyzing label ' num2str(i) ' of ' num2str(length(AND_mask)) ' (END)...']) + end + for j=i:length(AND_mask) + if i==j && ACh==2 + Fin_List{i}{j} = EDTI_Library.E_DTI_Net_An_dual_AND_analyze_ii(AND_mask([i j]),TInd, Lt, CLts); + else + Fin_List{i}{j} = EDTI_Library.E_DTI_Net_An_dual_AND_analyze(AND_mask([i j]),TInd, Lt, CLts); + end + end + end + end + + % From ExploreDTI: Helper function for connectivity analysis + function C = E_DTI_Get_IL_dist(mask) + CM = zeros(length(mask),3); + + L = length(mask); + + for i=1:L + + ind = find(mask{i}); + + [I,J,K] = ind2sub(size(mask{1}),ind); + + CM(i,1) = round(mean(I)); + CM(i,2) = round(mean(J)); + CM(i,3) = round(mean(K)); + + end + + C = repmat(nan,[L L]); + + for i=1:L + for j=i:L + + C(i,j) = sqrt(sum((CM(i,:)-CM(j,:)).^2)); + + end + end + + end + + % From ExploreDTI: Helper function for connectivity analysis + function FList = E_DTI_Net_An_dual_AND_analyze(AND_mask,TInd, Lt, CLts) + AND_IND = repmat(0>1,[length(TInd) length(AND_mask)]); + + for i=1:length(AND_mask) + AND_IND(:,i) = AND_mask{i}(TInd); + end + clear AND_mask TInd; + + FList = repmat(0>1,[Lt 1]); + + for i=1:Lt + FList(i) = all(any(AND_IND(CLts(i)+1:CLts(i+1),:),1),2); + end + clear AND_IND; + + FList_d = (1:Lt)'; + FList = FList_d(FList); + end + + % From ExploreDTI: Helper function for connectivity analysis + function FList = E_DTI_Net_An_dual_AND_analyze_ii(AND_mask,TInd, Lt, CLts) + AND_IND = repmat(0>1,[length(TInd) length(AND_mask)]); + + for i=1:length(AND_mask) + AND_IND(:,i) = AND_mask{i}(TInd); + end + clear AND_mask TInd; + + FList = repmat(0>1,[Lt 1]); + + for i=1:Lt + FList(i) = all(all(AND_IND(CLts(i)+1:CLts(i+1),:),1),2); + end + clear AND_IND; + + FList_d = (1:Lt)'; + FList = FList_d(FList); + end + + % From ExploreDTI: Helper function for connectivity analysis + function E_DTI_save_Network_analysis_temp_tracts(tract_f,F_List,out_f,L,sel_L,ACh) + [apf, anf, exte] = fileparts(tract_f); + if ACh==1 + fds = 'PASS'; + else + fds = 'END'; + end + + dir_temp = [out_f filesep 'temp_' fds '_' anf]; + + [suct, m1, m2] = mkdir(dir_temp); + + flag=0; + vlag=0; + + if suct==1 + + load(tract_f,'Tracts','VDims',... + 'TractMask','TractL','TractFE','TractFA','TractAng',... + 'TractGEO','TractMD','TractLambdas'); + + warning off all + load(tract_f,'TractAK','TractKA','TractRK','TractMK'); + warning on all + + if exist('TractAK','var') && exist('TractRK','var') ... + && exist('TractMK','var') && exist('TractKA','var') + flag=1; + TAK = TractAK; + TKA = TractKA; + TRK = TractRK; + TMK = TractMK; + clear TractAK TractKA TractRK TractMK; + end + + warning off all + load(tract_f,'TractsFOD'); + warning on all + + if exist('TractsFOD','var') + vlag=1; + TFOD = TractsFOD; + clear TractsFOD; + end + + + MDims = size(TractMask); + + for k=1:size(sel_L,1) + + i = sel_L(k,1); + j = sel_L(k,2); + + if i==j + nL = [L{i} '_tracts.mat']; + else + nL = [L{i} '_conn_' L{j} '_tracts.mat']; + end + + if i>j + i = sel_L(k,2); + j = sel_L(k,1); + end + + fname = [dir_temp filesep nL]; + + if ~isempty(F_List{i}{j}) + + % disp(['Saving fiber bundle ' num2str(k) ' of ' num2str(size(sel_L,1)) ' (' fds ') ...']) + + EDTI_Library.E_DTI_sve_NA_final(fname,Tracts(F_List{i}{j}),... + VDims,MDims,TractL(F_List{i}{j}),TractFE(F_List{i}{j}),... + TractFA(F_List{i}{j}),TractAng(F_List{i}{j}),... + TractGEO(F_List{i}{j}),TractMD(F_List{i}{j}),... + TractLambdas(F_List{i}{j})) + + if flag==1 + + TractAK = TAK(F_List{i}{j}); + TractRK = TRK(F_List{i}{j}); + TractMK = TMK(F_List{i}{j}); + TractKA = TKA(F_List{i}{j}); + + save(fname,'TractAK','TractKA','TractRK','TractMK','-append'); + + end + + if vlag==1 + + TractsFOD = TFOD(F_List{i}{j}); + save(fname,'TractsFOD','-append'); + + end + + + end + end + + end + end + + % From ExploreDTI: Helper function for connectivity analysis + function E_DTI_sve_NA_final(fname,Tracts,VDims,MDims,... + TractL,TractFE,TractFA,TractAng,TractGEO,TractMD,TractLambdas) + FList = (1:length(Tracts))'; + + TractMask = repmat(uint16(0),MDims); + + for i=1:length(FList) + + IND = unique(sub2ind(MDims,... + round(double(Tracts{FList(i)}(:,1))./VDims(1)),... + round(double(Tracts{FList(i)}(:,2))./VDims(2)),... + round(double(Tracts{FList(i)}(:,3))./VDims(3)))); + + TractMask(IND) = uint16(double(TractMask(IND)) + 1); + + end + + save(fname,'FList','Tracts','VDims','MDims',... + 'TractL','TractFE','TractFA','TractAng',... + 'TractGEO','TractMD','TractLambdas','TractMask') + end + + % From ExploreDTI: Helper function for connectivity analysis + function Out_names = E_DTI_get_NA_final_output_names(out_f,tract_f,Conn_Lab) + [apf, anf, exte] = fileparts(tract_f); + + for i=1:length(Conn_Lab) + Out_names{i} = [out_f filesep anf Conn_Lab{i} '.mat']; + end + end + + % From ExploreDTI: Helper function for connectivity analysis + function SA = E_DTI_Get_surface_estimate_of_mask(mask, Voxs) + sx = size(mask,1); + sy = size(mask,2); + sz = size(mask,3); + + Xdims = [1 sx]; + Ydims = [1 sy]; + Zdims = [1 sz]; + + x = single(Voxs(1):Voxs(1):Voxs(1)*Xdims(1,2)); + lx = length(x); + y = single(Voxs(2):Voxs(2):Voxs(2)*Ydims(1,2)); + ly = length(y); + z = single(Voxs(3):Voxs(3):Voxs(3)*Zdims(1,2)); + lz = length(z); + + ThreeD_X = single(repmat(x',[1 ly lz])); + ThreeD_Y = single(repmat(y,[lx 1 lz])); + p = single(zeros(1,1,lz)); + p(:,:,:) = z; + ThreeD_Z = single(repmat(p,[lx ly 1])); + + i_v = 0.5; + + [faces,verts] = isosurface(ThreeD_X,ThreeD_Y,ThreeD_Z,mask,i_v); + + a = verts(faces(:, 2), :) - verts(faces(:, 1), :); + b = verts(faces(:, 3), :) - verts(faces(:, 1), :); + c = cross(a, b, 2); + SA = 1/2 * sum(sqrt(sum(c.^2, 2))); + end + + % From ExploreDTI: Helper function for connectivity analysis + function M = E_DTI_NA_on_Vol(vol,T,VDims,MDims) + M = EDTI_Library.E_DTI_Get_Tracts_from_V(VDims,MDims,vol,T); + end + + % From ExploreDTI: Helper function for connectivity analysis + function TractVol = E_DTI_Get_Tracts_from_V(VDims,MDims,V,Tracts) + + V = double(V); + + res_method = 'cubic';% 'linear'; % or 'cubic' + + S = MDims; + + x = single(VDims(1):VDims(1):VDims(1)*S(1)); + lx = length(x); + y = single(VDims(2):VDims(2):VDims(2)*S(2)); + ly = length(y); + z = single(VDims(3):VDims(3):VDims(3)*S(3)); + lz = length(z); + + X = repmat(x',[1 ly lz]); + Y = repmat(y,[lx 1 lz]); + p = repmat(single(0),[1 1 lz]); + p(:,:,:) = z; + Z = repmat(p,[lx ly 1]); + clear x y z p lx ly lz; + + Lt = length(Tracts); + Lts = zeros(Lt,1); + + for i = 1:Lt + Lts(i,1)=size(Tracts{i},1); + end + CLts = [0 ;cumsum(Lts)]; + + Tmatx = repmat(single(0),[sum(Lts) 1]); + Tmaty = repmat(single(0),[sum(Lts) 1]); + Tmatz = repmat(single(0),[sum(Lts) 1]); + + for i=1:Lt + Tmatx(CLts(i)+1:CLts(i+1),1) = Tracts{i}(:,1); + Tmaty(CLts(i)+1:CLts(i+1),1) = Tracts{i}(:,2); + Tmatz(CLts(i)+1:CLts(i+1),1) = Tracts{i}(:,3); + end + + clear Tracts; + + Result = interpn(X,Y,Z,V,Tmatx,Tmaty,Tmatz,res_method); + + TractVol = cell(1,Lt); + for i=1:Lt + TractVol{i} = Result(CLts(i)+1:CLts(i+1),1); + end + end + + % From ExploreDTI: Helper function for connectivity analysis + function Vol = E_DTI_volume_NA(Tracts,VDims,MDims) + FList = (1:length(Tracts))'; + + TractMask = repmat(uint16(0),MDims); + + for i=1:length(FList) + + IND = unique(sub2ind(MDims,... + round(double(Tracts{FList(i)}(:,1))./VDims(1)),... + round(double(Tracts{FList(i)}(:,2))./VDims(2)),... + round(double(Tracts{FList(i)}(:,3))./VDims(3)))); + + TractMask(IND) = uint16(double(TractMask(IND)) + 1); + + end + + + Vol = TractMask>0; + Vol = sum(Vol(:))*prod(VDims); + end + + % From ExploreDTI: Helper function for connectivity analysis + function [Conn_matrices_1, CM1_ext] = E_DTI_Get_Network_analysis_conn_mat_1(tract_f,... + F_List,L,vol_f,dti_f, vol_ext,ACh,mask) + + if ACh==1 + suf_C = 'PASS'; + else + suf_C = 'END'; + end + + flag=0; + vlag=0; + + load(tract_f,'Tracts','VDims',... + 'TractMask','TractL','TractFE','TractFA','TractAng',... + 'TractGEO','TractMD','TractLambdas'); + + AT = length(Tracts); + + + warning off all + load(tract_f,'TractAK','TractKA','TractRK','TractMK'); + warning on all + + if exist('TractAK','var') && exist('TractRK','var') ... + && exist('TractMK','var') && exist('TractKA','var') + flag=1; + end + + warning off all + load(tract_f,'TractsFOD'); + warning on all + + if exist('TractsFOD','var') + vlag=1; + end + + Surf_A = zeros(length(mask),1); + + for i=1:length(mask) + if sum(mask{i}(:))>=1 + Surf_A(i,1) = EDTI_Library.E_DTI_Get_surface_estimate_of_mask(mask{i}, VDims); + else + Surf_A(i,1) = 1; + end + end + + MDims = size(TractMask); + + for i=1:15 + Conn_matrices_1{i} = repmat(nan,[length(L) length(L)]); + end + + CM1_ext{1} = ['_percentage_of_tracts_' suf_C]; + CM1_ext{2} = ['_average_tract_length_' suf_C]; + CM1_ext{3} = ['_tract_volume_' suf_C]; + CM1_ext{4} = ['_FA_' suf_C]; + CM1_ext{5} = ['_MD_' suf_C]; + CM1_ext{6} = ['_AD_' suf_C]; + CM1_ext{7} = ['_L2_' suf_C]; + CM1_ext{8} = ['_L3_' suf_C]; + CM1_ext{9} = ['_RD_' suf_C]; + CM1_ext{10} = ['_CL_' suf_C]; + CM1_ext{11} = ['_CP_' suf_C]; + CM1_ext{12} = ['_CS_' suf_C]; + CM1_ext{13} = ['_number_of_tracts_' suf_C]; + CM1_ext{14} = ['_density_weight_' suf_C]; + CM1_ext{15} = ['_binary_' suf_C]; + + + if flag==1 + + for i=1:4 + Conn_matrices_1{end+1} = repmat(nan,[length(L) length(L)]); + end + + CM1_ext{end+1} = ['_MK_' suf_C]; + CM1_ext{end+1} = ['_RK_' suf_C]; + CM1_ext{end+1} = ['_AK_' suf_C]; + CM1_ext{end+1} = ['_KA_' suf_C]; + + end + + if vlag==1 + + for i=1 + Conn_matrices_1{end+1} = repmat(nan,[length(L) length(L)]); + end + + CM1_ext{end+1} = ['_FOD_' suf_C]; + + end + + + Get_U_L = []; + + for i=1:length(L) + + % disp(['Calculating metrics for label ' num2str(i) ' of ' num2str(length(L)) ' (' suf_C ') ...']) + + for j=i:length(L) + + if ~isempty(F_List{i}{j}) + + Get_U_L = union(Get_U_L,F_List{i}{j}); + + Conn_matrices_1{1}(i,j) = length(F_List{i}{j}); + Conn_matrices_1{2}(i,j) = mean(cell2mat(TractL(F_List{i}{j}))); + Conn_matrices_1{3}(i,j) = EDTI_Library.E_DTI_volume_NA(Tracts(F_List{i}{j}),VDims,MDims); + Conn_matrices_1{4}(i,j) = mean(cell2mat(TractFA(F_List{i}{j})'))/sqrt(3); + Conn_matrices_1{5}(i,j) = mean(cell2mat(TractMD(F_List{i}{j})')); + eigv = cell2mat(TractLambdas(F_List{i}{j})'); + Conn_matrices_1{6}(i,j) = mean(eigv(:,1)); + Conn_matrices_1{7}(i,j) = mean(eigv(:,2)); + Conn_matrices_1{8}(i,j) = mean(eigv(:,3)); + Conn_matrices_1{9}(i,j) = (Conn_matrices_1{7}(i,j) + Conn_matrices_1{8}(i,j))/2; + GEO = cell2mat(TractGEO(F_List{i}{j})'); + Conn_matrices_1{10}(i,j) = mean(GEO(:,1)); + Conn_matrices_1{11}(i,j) = mean(GEO(:,2)); + Conn_matrices_1{12}(i,j) = mean(GEO(:,3)); + Conn_matrices_1{13}(i,j) = length(F_List{i}{j}); + Conn_matrices_1{14}(i,j) = (2/(Surf_A(i,1)+Surf_A(j,1)))*... + mean(1./cell2mat(TractL(F_List{i}{j}))); + Conn_matrices_1{15}(i,j) = 1; + + if flag==1 + + Conn_matrices_1{end-3}(i,j) = mean(cell2mat(TractMK(F_List{i}{j})')); + Conn_matrices_1{end-2}(i,j) = mean(cell2mat(TractRK(F_List{i}{j})')); + Conn_matrices_1{end-1}(i,j) = mean(cell2mat(TractAK(F_List{i}{j})')); + Conn_matrices_1{end}(i,j) = mean(cell2mat(TractKA(F_List{i}{j})')); + + end + + if vlag==1 + + Conn_matrices_1{end}(i,j) = mean(cell2mat(TractsFOD(F_List{i}{j})')); + + end + + end + end + end + + Conn_matrices_1{1} = Conn_matrices_1{1}/length(Get_U_L); + Conn_matrices_1{15} = ~isnan(Conn_matrices_1{15}); + + if ~isempty(vol_f) + + flag_ = 1; + + if exist(vol_f,'file')~=2 + disp('Warning: it seems that the following file cannot be read:') + disp(vol_f) + flag_ = 0; + end + + try + [vol, VDims_vol] = EDTI_Library.E_DTI_read_nifti_file(vol_f); + catch me + disp('Error:') + disp(me.message) + flag_ = 0; + end + + if ndims(vol)<3 || ndims(vol)>4 + flag_ = 0; + end + + load(dti_f,'MDims','VDims') + + if ~all(MDims==size(vol(:,:,:,1))) + flag_ = 0; + end + + if ~all(VDims==VDims_vol) + flag_ = 0; + end + + if flag_==1 + + vol = double(vol); + + for v=1:size(vol,4) + + % disp(['Calculating ''volume'' metrics for volume ' num2str(v) ' of ' num2str(size(vol,4)) ' (' suf_C ') ...']) + + Conn_matrices_1{end+1} = repmat(nan,[length(L) length(L)]); + + if size(vol,4)==1 + CM1_ext{end+1} = [vol_ext(1:end-4) '_' suf_C]; + else + CM1_ext{end+1} = [vol_ext(1:end-4) '_' num2str(v) '_' suf_C]; + end + + for i=1:length(L) + for j=i:length(L) + + if ~isempty(F_List{i}{j}) + + TV = EDTI_Library.E_DTI_NA_on_Vol(vol(:,:,:,v),... + Tracts(F_List{i}{j}),VDims,MDims); + + Conn_matrices_1{end}(i,j) = ... + mean(cell2mat(TV(1:end)')); + + end + end + end + + end + end + end + + for p=1:length(Conn_matrices_1) + for i=1:size(Conn_matrices_1{p},1) + for j=i:size(Conn_matrices_1{p},2) + + Conn_matrices_1{p}(j,i) = Conn_matrices_1{p}(i,j); + + end + end + end + end + % This function is used to track multiple FODs (mFOD). Written by A. De % Luca function WholeBrainTracking_mDRL_fast_exe(f_in, suffix, f_out, p, use_linear_or_majority) @@ -8549,7 +9427,7 @@ function WholeBrainTracking_mDRL_fast_exe(f_in, suffix, f_out, p, use_linear_or_ end - % Wraooer around the Gibbs ringing correction + % Wrapper around the Gibbs ringing correction function GibbsRingingCorrection(fnp_in,fnp_out) try diff --git a/ExploreDTIInterface/mk_init.m b/ExploreDTIInterface/mk_init.m index 23e1578..196a2e1 100644 --- a/ExploreDTIInterface/mk_init.m +++ b/ExploreDTIInterface/mk_init.m @@ -1,13 +1 @@ -%%%$ Included in MRIToolkit (https://github.com/delucaal/MRIToolkit) %%% %%% Alberto De Luca - alberto@isi.uu.nl $%%% %%% Distributed under the terms of LGPLv3 %%% - - - -% A. De Luca -% Sub-toolbox specific initialization -% 09/11/2019: creation - v1.0 -global MRIToolkit; -MRIToolkit.exploredtiinterface_version = 1.0; - -addpath(get_executed_file_path()) -addpath(fullfile(get_executed_file_path(),'Auxiliaries')) -addpath(fullfile(get_executed_file_path(),'Auxiliaries','MEX')) \ No newline at end of file +%%%$ Included in MRIToolkit (https://github.com/delucaal/MRIToolkit) %%% %%% Alberto De Luca - alberto@isi.uu.nl $%%% %%% Distributed under the terms of LGPLv3 %%% % A. De Luca % Sub-toolbox specific initialization % 09/11/2019: creation - v1.0 % 18/12/2024: creation - v1.0 global MRIToolkit; MRIToolkit.exploredtiinterface_version = 1.1; addpath(get_executed_file_path()) addpath(fullfile(get_executed_file_path(),'Auxiliaries')) addpath(fullfile(get_executed_file_path(),'Auxiliaries','MEX')) \ No newline at end of file diff --git a/MRIToolkitInit.m b/MRIToolkitInit.m index a0f1a1b..d2d8912 100644 --- a/MRIToolkitInit.m +++ b/MRIToolkitInit.m @@ -13,6 +13,7 @@ % 05/05/2020: Added suppot for automatic version control % 13/08/2021: v1.5 % 12/01/2024: v1.6 +% 18/12/2024: v1.7 function MRIToolkitInit(SelectedToolboxes) run_folder = mfilename('fullpath'); if(~isempty(run_folder)) @@ -114,8 +115,8 @@ function MRIToolkitInit(SelectedToolboxes) git_version = fgetl(git_commit); fclose(git_commit); elseif(isdeployed) - git_version = 'v1.6-cmdline'; + git_version = 'v1.7-cmdline'; else - git_version = 'v1.6-nongit'; + git_version = 'v1.7-nongit'; end end diff --git a/Neuroimaging/Neuro.m b/Neuroimaging/Neuro.m index ba7a824..1265489 100644 --- a/Neuroimaging/Neuro.m +++ b/Neuroimaging/Neuro.m @@ -187,6 +187,38 @@ function BrainExtractionWithRegistration(varargin) MRTQuant.WriteNifti(in_file,output); rmdir(temp_direc,'s'); end + + % Reconstruct structural connectivity given whole brain tractography + % and corresponding brain parcellations in native space. + function StructuralConnectivityAnalysis(varargin) + coptions = varargin; + mat_file = GiveValueForName(coptions,'mat_file'); + if(isempty(mat_file)) + error('Missing mandatory argument mat_file'); + end + tract_file = GiveValueForName(coptions,'tract_file'); + if(isempty(tract_file)) + error('Missing mandatory argument tract_file'); + end + label_file = GiveValueForName(coptions,'label_file'); + if(isempty(label_file)) + error('Missing mandatory argument label_file'); + end + labels_txt = GiveValueForName(coptions,'labels_txt'); + if(isempty(labels_txt)) + error('Missing mandatory argument labels_txt'); + end + output = GiveValueForName(coptions,'output'); + if(isempty(output)) + error('Missing mandatory argument output'); + end + if(exist(output,'dir') < 1) + mkdir(output); + end + EDTI_Library.E_DTI_Network_analysis_from_ROI_L({mat_file},{tract_file},{label_file},labels_txt,1,... + 'pass',output,0); + + end end end diff --git a/Neuroimaging/mk_init.m b/Neuroimaging/mk_init.m index e9392a0..58488b7 100644 --- a/Neuroimaging/mk_init.m +++ b/Neuroimaging/mk_init.m @@ -1 +1 @@ -%%%$ Included in MRIToolkit (https://github.com/delucaal/MRIToolkit) %%% %%% Alberto De Luca - alberto@isi.uu.nl $%%% %%% Distributed under the terms of LGPLv3 %%% % A. De Luca % Sub-toolbox specific initialization % 20/10/2020: creation - v1.1 global MRIToolkit; MRIToolkit.neuro_version = 1.1; addpath(get_executed_file_path()) \ No newline at end of file +%%%$ Included in MRIToolkit (https://github.com/delucaal/MRIToolkit) %%% %%% Alberto De Luca - alberto@isi.uu.nl $%%% %%% Distributed under the terms of LGPLv3 %%% % A. De Luca % Sub-toolbox specific initialization % 20/10/2020: creation - v1.1 % 18/12/2024: v1.2 global MRIToolkit; MRIToolkit.neuro_version = 1.2; addpath(get_executed_file_path()) \ No newline at end of file diff --git a/README.md b/README.md index 69d31e1..dad8442 100644 --- a/README.md +++ b/README.md @@ -4,67 +4,12 @@

-# [MRIToolkit](https://github.com/delucaal/MRIToolkit) [Update 11-12-2024] -New in this version [1.6b]: -- Bug fixes -- Updated compatibility to the latest version of WhiteMatterAnalysis -- Initial support for tractography with Scilpy tools +# [MRIToolkit](https://github.com/delucaal/MRIToolkit) [Update 11-12-2024 v1.7] -Update 17-02-2023: -- This is the last stable release before a major change in structure. From the next release, all main functionalities will be come available as standalone commands and be designed to fully rely on Niftis only. -- Fixes to the handling of Nifti headers -- Automated brain extraction using image registration -- added new class MRT_Library that will collect all own methods currently in MRTTrack and MRTQuant (transition in progress). This will clean MRTTrack and MRTQuant which are meant as documented interfaces. -- Dramatic speed up of GRL and mFOD (major revision of MRTTrack.PerformDeconv) -- Changes to MRTQuant.LoadNifti and SaveNifti, which now allows to preserve the original header of NIFTIs (if required). In future releases, no modifications of the headers will be required. -- Update of Neuro.m to support the latest version of Elastix -- Compressed niftis are now (un)compressed in a temporary folder -- Added functions to track multiple FODs simultaneously -- Changes to the automatic handling of Q-S form of NIFTIs -- New mFOD fiber tracking -- New TerminateTractsWithFraction approach with 1 voxel tolerance -- Moved Trackers to aanother location -- Fixes to SHPrecomp -- Fixes to MRIToolkitInit -- Tractography can now use parallel computing -- Handling of NIFTI scaling factors -- Failsafe mode for registration -- Support for FSL Eddy -- Support for FOD normalization -- One call command to run GRL -- Fix for MSCSD (probably still needs amendment on the detection of data shells) -- Fix fot CAT12 processing -- command line fixes -- standardized preproc includes conformation -- fixes to conformation -- support for offset in DKI for GRL -- fixes to DKI export - - -Update 09-2021: -- [NEW] Checkout the [tutorial](https://github.com/delucaal/MRIToolkit/tree/master/GettingStarted/Diffusion/7_WMA) on how to install/use automated tractography clustering (WMA) algorithm with CSD/GRL/mFOD -- Fixes to MSCSD-related functions -- Fixes to WM automated clustering -- Support for Q-form in MRTTrack.ConformSpatialDimensions() -- New FOD scaling mode for GRL/mFOD -- Renaming of some command line functions -- Added an option to customise the number of dRL iterations -- New functions to estimate SNR from b=0s/mm2 images (MRTQuant) - - -Update 17-02-2021: -- Run the CAT12 automatic pipeline for T1 images -- Code re-organization into two main classes: MRTQuant (preprocessing/DTI/DKI) and MRTrack (Tractography related) -- Support for automatic fiber clustering (see below for reference) -- Early support for integration with Python (needed for the point above) -- Support for VTK poly data 4.2 and 5.1 (also needed for the clustering) -- Added a robust option to GRL/mFOD deconvolution -- Initial support for storing the NIFTI Q/S form (to improve interoperability with other tools, not implemented yet) -- Integration with CAT12 ## What is it? -MRIToolkit is a set of [command line tools](https://github.com/delucaal/MRIToolkit/tree/master/CommandLine) and a MATLAB (R) toolbox/library to process (diffusion) magnetic resonance imaging (MRI) data. Binaries of the command line version will be provided soon! +MRIToolkit is a set of [command line tools](https://github.com/delucaal/MRIToolkit/tree/master/CommandLine) and a MATLAB (R) toolbox to process (diffusion) magnetic resonance imaging (MRI) data. Binaries of the command line version will be provided at some point! -The idea behind MRIToolkit is to integrate [my research output](https://www.umcutrecht.nl/en/research/researchers/de-luca-alberto-a) with existing state-of-the-art methods for (diffusion) MRI processing. +The idea behind MRIToolkit is to integrate [my research output](https://research.umcutrecht.nl/researchers/de-luca/) with existing state-of-the-art methods for (diffusion) MRI processing. ## Where do I find it? - The MATLAB (R) toolbox is available [here](https://github.com/delucaal/MRIToolkit) on Github!