123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264 |
- % thePlotSpm
- %
- % This is a Matlab script that uses custom code and spm12 routines to plot
- % a version of THE PLOT from Jonathan Power. See:
- % - https://www.jonathanpower.net/2017-ni-the-plot.html
- % - https://www.ncbi.nlm.nih.gov/pubmed/27510328
- % (code has been partly adopted from Power's script)
- %
- % THE PLOT is a useful visual quality checking tool for fMRI timeseries
- % data. It is a 2D plot of voxel intensities over time, with voxels ordered
- % in bins (GM, WM, CSF) to indicate some directionality in the data. Motion
- % and physiological noise are more easily visually inspected with this
- % tool, especially if viewed next to other quality traces, like FD, DVARS
- % or physiological recordings.
- % My script does the following:
- % - Run preprocessing on input data:
- % - Checks if T1w segments already exist.
- % - If not, realigns all images in time series to first functional
- % image, coregisters (estimate) T1w image to first functional image,
- % segments the coregistered T1w image into GM, WM and CSF
- % compartments, reslices all segments into functional image grid,
- % smooths realigned functional images, smooths unprocessed functional
- % images.
- % - Calculates framewise displacement
- % - Creates brain mask (GM+WM+CSF)
- % - Calculates BOLD percentage signal changes for masked brain voxels
- % - Creates figure with FD and THEPLOT (with blue and red lines
- % separating GM (top), WM (middle), and CSF(bottom))
- %
- % INPUTS:
- % structural_fn - filename of T1-weighted structural scan, a 4D nifti
- % file is expected (e.g. 'mprage.nii')
- % funcional4D_fn - filename of main functional scan, a 4D nifti file is
- % expected (e.g. 'rest.nii')
- %
- % DEPENDENCIES:
- % thePlotSpmPreproc.m - preprocesses data (functional realignment, coregistration of
- % T1w to first functional image, segmentation into
- % tissue types, reslicing to functional grid)
- % createBinarySegments.m - constructs 3D binary images for GM, WM and CSF based
- % on the relative value of GM/WM/CSF tissue probability
- % maps per voxel.
- %
- % NOTES/ASSUMPTIONS:
- % - This script makes use of spm12 batch scripting
- % - This script assumes that image files are in the Nifti file format
- % - This script was developed in Matlab 2016b on a Mac OS Sierra, and not
- % tested on other platforms (it wasn't even tested much on my Mac)
- % - This is a simplified, Matlab-SPM-only version of Power's THEPLOT and
- % excludes things like erosion of masks and segmentation into more bins
- % (done in freesurfer).
- % - I might have made some mistakes in the code, or it could be improved in
- % many ways. Please contribute!
- %
- % 2018-08-17 Remi Gau
- % This seems to work on windows 10 with matlab 2017a.
- % While looking around I realized that Cyril Pernet has a function
- % (spmup_timeseriesplot) to plot this (and do quite a few other things as
- % well) in his spmup toolbox: % https://github.com/CPernet/spmup
- % But I like your script a lot as well as it is more stand-alone and
- % therefore more pedagogical.
- %
- %__________________________________________________________________________
- % Copyright (C) 2018 - Stephan Heunis
- clear
- clc
- %% User defined variables:
- % -------------------------------------------------------------------------
- data_dir = ''; % e.g. '/users/me/matlab/data/subj1'
- functional4D_fn = [data_dir filesep '']; % e.g. [data_dir filesep 'rest.nii']
- structural_fn = [data_dir filesep '']; % e.g. [data_dir filesep 'mprage.nii']
- spm_dir = ''; % e.g. '/users/me/matlab/spm12'
- fwhm = 6; % for preprocessing smoothing steps
- use_processed = 0; % use realigned data for the plot? yes = 1
- intensity_scale = [-6 6]; % scaling for plot image intensity, see what works
- davg_def = 50; %Average distance to the surface of the brain ; 50mm = assumed brain radius; from article
- % -------------------------------------------------------------------------
- % RG
- % -------------------------------------------------------------------------
- fwhm = 6;
- [spm_dir, ~, ~] = fileparts(which('spm'));
- data_dir = 'D:';
- subj = 'sub-01';
- % structural image
- structural_fn = spm_select('FPList', fullfile(data_dir , subj , 'anat'), ['^' subj '_T1w.*nii$']);
- % 4D BOLD time series
- functional4D_fn = spm_select('FPList', fullfile(data_dir , subj , 'func'), ['^' subj '_task-.*_bold*.nii$']);
- use_processed = 0; % use realigned data for the plot? yes = 1
- intensity_scale = [-6 6]; % scaling for plot image intensity, see what works
- davg_def = 65;
- % -------------------------------------------------------------------------
- %% Get image information
- func_spm = spm_vol(functional4D_fn);
- tsize = size(func_spm);
- Nt = tsize(1);
- Ni= func_spm(1).dim(1);
- Nj= func_spm(1).dim(2);
- Nk= func_spm(1).dim(3);
- %% Data check
- % Preprocess structural and functional data. First check if a resliced grey
- % matter image is present in the specified data directory. If true, assume
- % that all preprocessing has been completed by thePlotSpmPreproc, declare
- % appropriate variables and
- [dir, file, ext] = fileparts(structural_fn);
- if exist([dir filesep 'rc1' file ext], 'file')
- disp('Preproc done!')
- preproc_data = struct;
- preproc_data.rstructural_fn = fullfile(dir , ['r' file ext]);
- preproc_data.rgm_fn = fullfile(dir , ['rc1' file ext]);
- preproc_data.rwm_fn = fullfile(dir , ['rc2' file ext]);
- preproc_data.rcsf_fn = fullfile(dir , ['rc3' file ext]);
-
- [dir, file, ext] = fileparts(functional4D_fn);
- preproc_data.rfunctional_fn = fullfile(dir , ['r' file ext]);
- preproc_data.srfunctional_fn = fullfile(dir , ['sr' file ext]);
- preproc_data.sfunctional_fn = fullfile(dir , ['s' file ext]);
- preproc_data.mp_fn = fullfile(dir , ['rp_' file '.txt']);
- preproc_data.MP = load(preproc_data.mp_fn);
- else
- preproc_data = thePlotSpmPreproc(functional4D_fn, structural_fn, fwhm, spm_dir);
- end
- % Additional check to make sure that the functional images and the resliced
- % TPMs have same dimensions
- files = {};
- files{1,1} = preproc_data.rgm_fn;
- files{2,1} = [preproc_data.rfunctional_fn ',1'];
- files = char(files);
- spm_check_orientations(spm_vol(files));
- %% Try to compute the average surface to the brain
- % from motion finger print functions and script (mw_mfp.m) from Marko Wilke
- % https://www.medizin.uni-tuebingen.de/kinder/en/research/neuroimaging/software/
- % see http://www.dx.doi.org/10.1016/j.neuroimage.2011.10.043
- % and http://www.dx.doi.org/10.1371/journal.pone.0106498
- [dir, file, ext] = fileparts(structural_fn);
- try
- if isempty(spm_select('FPList',dir,'^rc1.*\.surf\.gii$'))
- % create a surface of the brain using the TPM if we don't have one
- files = spm_select('FPList',dir,'^rc[12].*\.nii$');
- spm_surf(files,2);
- clear files
- end
- % compute the average surface to the brain
- FV = gifti(spm_select('FPList',dir,'^rc1.*\.surf\.gii$'));
- center = FV.vertices(FV.faces(:,:),:);
- center = reshape(center, [size(FV.faces,1) 3 3]);
- center = squeeze(mean(center,2));
- ori_dist = sqrt(sum((center.*-1).^2,2))';
- davg = mean(ori_dist);
- clear ori_dist center FV
- catch
- davg = davg_def;
- end
- fprintf(' Average distance to the cortex surface: %0.2f ', davg)
- %% Calculate Framewise displacement:
- % - First demean and detrend the movement parameters and
- % - then calculate FD.
- % See also spmup_FD from https://github.com/CPernet/spmup
- MP2 = preproc_data.MP - repmat(mean(preproc_data.MP, 1), Nt,1);
- X = (1:Nt)';
- X2 = X - mean(X);
- X3 = [ones(Nt,1) X2];
- b = X3\MP2;
- MP_corrected = MP2 - X3(:, 2)*b(2,:);
- MP_mm = MP_corrected;
- MP_mm(:,4:6) = MP_mm(:,4:6)*davg;
- MP_diff = [zeros(1, 6); diff(MP_mm)];
- FD = sum(abs(MP_diff),2);
- %% Calculate BOLD percentage signal change (PSC)
- % - First create binary masks for GM, WM and CSF,
- % - Then combine them to get a brain mask for which calculations are run
- % - Then compute PSC for unprocessed or realigned data (user selection)
- % - Then plot everything
- [GM_img_bin, WM_img_bin, CSF_img_bin] = createBinarySegments(preproc_data.rgm_fn, preproc_data.rwm_fn, preproc_data.rcsf_fn, 0.1);
- I_GM = find(GM_img_bin);
- I_WM = find(WM_img_bin);
- I_CSF = find(CSF_img_bin);
- mask_reshaped = GM_img_bin | WM_img_bin | CSF_img_bin;
- I_mask = find(mask_reshaped);
- Nmaskvox = numel(I_mask);
- %% Load and detrend data
- % Use realigned or unprocessed data for plot, user-selected
- if use_processed == 1
- F_img = spm_read_vols(spm_vol(preproc_data.srfunctional_fn));
- else
- F_img = spm_read_vols(spm_vol(preproc_data.sfunctional_fn));
- end
- % Remove linear and polynomial trends from data
- F_2D = reshape(F_img, Ni*Nj*Nk, Nt);
- X_design = [ (1:Nt)' ((1:Nt).^2/(Nt^2))' ((1:Nt).^3/(Nt^3))'];
- X_design = X_design - mean(X_design);
- X_design = [X_design ones(Nt,1)];
- betas = X_design\F_2D';
- F_detrended = F_2D' - X_design(:, 1:(end-1))*betas(1:(end-1), :);
- F_detrended = F_detrended';
- % Mask, mean and PSC
- F_masked = F_detrended(I_mask, :);
- F_mean = mean(F_masked, 2);
- F_masked_psc = 100*(F_masked./repmat(F_mean, 1, Nt)) - 100;
- F_masked_psc(isnan(F_masked_psc))=0;
- F_psc_img = zeros(Ni, Nj, Nk, Nt);
- F_2D_psc = reshape(F_psc_img, Ni*Nj*Nk, Nt);
- F_2D_psc(I_mask, :) = F_masked_psc;
- F_psc_img = reshape(F_2D_psc, Ni, Nj, Nk, Nt);
- %% Figure
- % Create image to plot by concatenation
- GM_img = F_2D_psc(I_GM, :);
- WM_img = F_2D_psc(I_WM, :);
- CSF_img = F_2D_psc(I_CSF, :);
- all_img = [GM_img; WM_img; CSF_img];
- % Identif limits between the different tissue compartments
- line1_pos = numel(I_GM);
- line2_pos = numel(I_GM) + numel(I_WM);
- figure
- fontsizeL = 17;
- fontsizeM = 15;
- % plot functional data
- ax1 = subplot(5,1,[2:5]);
- imagesc(ax1, all_img);
- colormap(gray);
- caxis(intensity_scale);
- title(ax1, 'thePlotSpm','fontsize',fontsizeL)
- ylabel(ax1, 'Voxels','fontsize',fontsizeM)
- xlabel(ax1, 'fMRI volumes','fontsize',fontsizeM)
- hold on;
- line([1 Nt],[line1_pos line1_pos], 'Color', 'b', 'LineWidth', 2 )
- line([1 Nt],[line2_pos line2_pos], 'Color', 'r', 'LineWidth', 2 )
- hold off;
- % plot Framewise displacement
- ax2 = subplot(5,1,1);
- plot(ax2, FD, 'LineWidth', 2); grid;
- axis tight
- set(ax2,'Xticklabel',[]);
- title(ax2, 'FD','fontsize',fontsizeL)
- ylabel(ax2, 'mm','fontsize',fontsizeM)
|