Skip to content

[ENH] add resize_image function to library #297

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Feb 17, 2021
Merged

[ENH] add resize_image function to library #297

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
58 changes: 29 additions & 29 deletions lib/utils/plot_power_spectra_of_GLM_residuals.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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;

Expand All @@ -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');
Expand All @@ -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
Expand Down Expand Up @@ -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);
Expand All @@ -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');
Expand Down
149 changes: 149 additions & 0 deletions lib/utils/resize_img.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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];