diff --git a/lib/utils/plot_power_spectra_of_GLM_residuals.m b/lib/utils/plot_power_spectra_of_GLM_residuals.m index b99acdec..e4b55182 100644 --- a/lib/utils/plot_power_spectra_of_GLM_residuals.m +++ b/lib/utils/plot_power_spectra_of_GLM_residuals.m @@ -2,21 +2,21 @@ function plot_power_spectra_of_GLM_residuals(path_to_results, TR, cutoff_freq, a % -By Wiktor Olszowy, University of Cambridge, olszowyw@gmail.com % - % -Written following study + % -Written following study % 'Accurate autocorrelation modeling substantially improves fMRI reliability' % -https://www.nature.com/articles/s41467-019-09230-w.pdf % -May 2018 % - % -Given fMRI task results in AFNI, FSL or SPM, + % -Given fMRI task results in AFNI, FSL or SPM, % this script plots power spectra of GLM residuals. - % -If there is strong structure visible in the GLM residuals + % -If there is strong structure visible in the GLM residuals % the power spectra are not flat), the first level results are likely confounded. % -tested on Linux % -you need on your path >= MATLAB 2017b, AFNI and FSL % -specify the default values for the cutoff frequency used by the high-pass filter, - % -for the assumed experimental design frequency + % -for the assumed experimental design frequency % and for the true experimental design frequency; % -10 chosen, as it is beyond the plotted frequencies @@ -29,11 +29,11 @@ function plot_power_spectra_of_GLM_residuals(path_to_results, TR, cutoff_freq, a if ~exist('true_exper_freq', 'var') true_exper_freq = 10; end - + Fontsize = 10; - % -Fast Fourier Transform (FFT) will pad the voxel-wise time series - % to that length with trailing zeros (if no. of time points lower) + % -Fast Fourier Transform (FFT) will pad the voxel-wise time series + % to that length with trailing zeros (if no. of time points lower) % or truncate to that length (if no. of time points higher) fft_n = 512; @@ -49,13 +49,13 @@ function plot_power_spectra_of_GLM_residuals(path_to_results, TR, cutoff_freq, a system(['3dcalc -a ' AFNI_res4d_name ' -expr "a" -prefix res4d.nii']); res4d = niftiread('res4d.nii'); - % -for FSL + % -for FSL elseif exist('stats/res4d.nii.gz', 'file') == 2 res4d = niftiread('stats/res4d.nii.gz'); elseif exist('stats/res4d.nii', 'file') == 2 res4d = niftiread('stats/res4d.nii'); - % -for SPM + % -for SPM elseif exist('Res_0001.nii', 'file') == 2 SPM_res4d_name = dir('Res_*.nii'); @@ -73,15 +73,15 @@ function plot_power_spectra_of_GLM_residuals(path_to_results, TR, cutoff_freq, a end system(['fslmerge -t res4d ' SPM_res4d_all]); res4d = niftiread('res4d.nii.gz'); - + end else disp(['No GLM residuals found! ', ... - 'If you run SPM, remember to put command ', ... - 'VRes = spm_write_residuals(SPM, NaN) at the end of the SPM script. ', ... - 'Otherwise, SPM by default deletes the GLM residuals.']); + 'If you run SPM, remember to put command ', ... + 'VRes = spm_write_residuals(SPM, NaN) at the end of the SPM script. ', ... + 'Otherwise, SPM by default deletes the GLM residuals.']); return end @@ -134,7 +134,7 @@ function plot_power_spectra_of_GLM_residuals(path_to_results, TR, cutoff_freq, a % -make the plot figure('rend', 'painters', 'pos', [0 0 600 400], 'Visible', 'on'); - hold on + hold on; f = linspace(0, 0.5 / TR, 257); max_y = max(power_spectra_of_GLM_residuals); @@ -160,28 +160,28 @@ function plot_power_spectra_of_GLM_residuals(path_to_results, TR, cutoff_freq, a location = 'southeast'; if (power_spectra_of_GLM_residuals(257) < 0) || (max_y > 2) - location = 'northeast'; + location = 'northeast'; end - + to_plot = [h1 h10 h4 h6 h8]; - legend_content = {... - 'Actual power spectrum', ... - 'Ideal power spectrum', ... - 'High pass filter frequency cutoff', ... - 'Assumed design frequency', ... - 'True design frequency'}; - + legend_content = { ... + 'Actual power spectrum', ... + 'Ideal power spectrum', ... + 'High pass filter frequency cutoff', ... + 'Assumed design frequency', ... + 'True design frequency'}; + if any([cutoff_freq, assumed_exper_freq, true_exper_freq] == 10) to_plot = [h1 h10]; - legend_content = {... - 'Actual power spectrum', ... - 'Ideal power spectrum'}; + legend_content = { ... + 'Actual power spectrum', ... + 'Ideal power spectrum'}; end legend(to_plot, legend_content, ... - 'box', 'off', ... - 'FontSize', Fontsize, ... - 'Location', location); + 'box', 'off', ... + 'FontSize', Fontsize, ... + 'Location', location); figname = 'power_spectra_of_GLM_residuals'; print (figname, '-dpng'); diff --git a/lib/utils/resize_img.m b/lib/utils/resize_img.m new file mode 100644 index 00000000..fdd5407d --- /dev/null +++ b/lib/utils/resize_img.m @@ -0,0 +1,149 @@ +function resize_img(imnames, Voxdim, BB, ismask) + % resize_img -- resample images to have specified voxel dims and BBox + % + % resize_img(imnames, voxdim, bb, ismask) + % + % Output images will be prefixed with 'r', and will have voxel dimensions + % equal to voxdim. Use NaNs to determine voxdims from transformation matrix + % of input image(s). + % If bb == nan(2,3), bounding box will include entire original image + % Origin will move appropriately. Use world_bb to compute bounding box from + % a different image. + % + % Pass ismask=true to re-round binary mask values (avoid + % growing/shrinking masks due to linear interp) + % + % See also voxdim, world_bb + + % Based on John Ashburner's reorient.m + % http://www.sph.umich.edu/~nichols/JohnsGems.html#Gem7 + % http://www.sph.umich.edu/~nichols/JohnsGems5.html#Gem2 + % Adapted by Ged Ridgway -- email bugs to drc.spm@gmail.com + + % This version doesn't check spm_flip_analyze_images -- the handedness of + % the output image and matrix should match those of the input. + + % DONWLOADED FROM ON THE 2021 02 17 + % https://blogs.warwick.ac.uk/nichols/entry/spm5_gem_3/ + + % Check spm version: + if exist('spm_select', 'file') % should be true for spm5 + spm5 = 1; + elseif exist('spm_get', 'file') % should be true for spm2 + spm5 = 0; + else + error('Can''t find spm_get or spm_select; please add SPM to path'); + end + + % spm_defaults; + + % prompt for missing arguments + if ~exist('imnames', 'var') || isempty(char(imnames)) + if spm5 + imnames = spm_select(inf, 'image', 'Choose images to resize'); + else + imnames = spm_get(inf, 'img', 'Choose images to resize'); + end + end + % check if inter fig already open, don't close later if so... + Fint = spm_figure('FindWin', 'Interactive'); + Fnew = []; + if ~exist('Voxdim', 'var') || isempty(Voxdim) + Fnew = spm_figure('GetWin', 'Interactive'); + Voxdim = spm_input('Vox Dims (NaN for "as input")? ', ... + '+1', 'e', '[nan nan nan]', 3); + end + if ~exist('BB', 'var') || isempty(BB) + Fnew = spm_figure('GetWin', 'Interactive'); + BB = spm_input('Bound Box (NaN => original)? ', ... + '+1', 'e', '[nan nan nan; nan nan nan]', [2 3]); + end + if ~exist('ismask', 'var') + ismask = false; + end + if isempty(ismask) + ismask = false; + end + + % reslice images one-by-one + vols = spm_vol(imnames); + for V = vols' + % (copy to allow defaulting of NaNs differently for each volume) + voxdim = Voxdim; + bb = BB; + % default voxdim to current volume's voxdim, (from mat parameters) + if any(isnan(voxdim)) + vprm = spm_imatrix(V.mat); + vvoxdim = vprm(7:9); + voxdim(isnan(voxdim)) = vvoxdim(isnan(voxdim)); + end + voxdim = voxdim(:)'; + + mn = bb(1, :); + mx = bb(2, :); + % default BB to current volume's + if any(isnan(bb(:))) + vbb = world_bb(V); + vmn = vbb(1, :); + vmx = vbb(2, :); + mn(isnan(mn)) = vmn(isnan(mn)); + mx(isnan(mx)) = vmx(isnan(mx)); + end + + % voxel [1 1 1] of output should map to BB mn + % (the combination of matrices below first maps [1 1 1] to [0 0 0]) + mat = spm_matrix([mn 0 0 0 voxdim]) * spm_matrix([-1 -1 -1]); + % voxel-coords of BB mx gives number of voxels required + % (round up if more than a tenth of a voxel over) + imgdim = ceil(mat \ [mx 1]' - 0.1)'; + + % output image + VO = V; + [pth, nam, ext] = fileparts(V.fname); + VO.fname = fullfile(pth, ['r' nam ext]); + VO.dim(1:3) = imgdim(1:3); + VO.mat = mat; + VO = spm_create_vol(VO); + spm_progress_bar('Init', imgdim(3), 'reslicing...', 'planes completed'); + for i = 1:imgdim(3) + M = inv(spm_matrix([0 0 -i]) * inv(VO.mat) * V.mat); + img = spm_slice_vol(V, M, imgdim(1:2), 1); % (linear interp) + if ismask + img = round(img); + end + spm_write_plane(VO, img, i); + spm_progress_bar('Set', i); + end + spm_progress_bar('Clear'); + end + % call spm_close_vol if spm2 + if ~spm5 + spm_close_vol(VO); + end + if isempty(Fint) && ~isempty(Fnew) + % interactive figure was opened by this script, so close it again. + close(Fnew); + end + disp('Done.'); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function bb = world_bb(V) + % world-bb -- get bounding box in world (mm) coordinates + + d = V.dim(1:3); + % corners in voxel-space + c = [1 1 1 1 + 1 1 d(3) 1 + 1 d(2) 1 1 + 1 d(2) d(3) 1 + d(1) 1 1 1 + d(1) 1 d(3) 1 + d(1) d(2) 1 1 + d(1) d(2) d(3) 1]'; + % corners in world-space + tc = V.mat(1:3, 1:4) * c; + + % bounding box (world) min and max + mn = min(tc, [], 2)'; + mx = max(tc, [], 2)'; + bb = [mn; mx];