diff --git a/dependencies/MNIBrainMask.mat b/dependencies/MNIBrainMask.mat new file mode 100644 index 0000000..eba3a58 Binary files /dev/null and b/dependencies/MNIBrainMask.mat differ diff --git a/dependencies/VOI_PCC_1.mat b/dependencies/VOI_PCC_1.mat new file mode 100644 index 0000000..a80c0c9 Binary files /dev/null and b/dependencies/VOI_PCC_1.mat differ diff --git a/dependencies/cbrewer/cbrewer/cbrewer.m b/dependencies/cbrewer/cbrewer/cbrewer.m new file mode 100644 index 0000000..26be891 --- /dev/null +++ b/dependencies/cbrewer/cbrewer/cbrewer.m @@ -0,0 +1,128 @@ +function [colormap]=cbrewer(ctype, cname, ncol, interp_method) +% +% CBREWER - This function produces a colorbrewer table (rgb data) for a +% given type, name and number of colors of the colorbrewer tables. +% For more information on 'colorbrewer', please visit +% http://colorbrewer2.org/ +% +% The tables were generated from an MS-Excel file provided on the website +% http://www.personal.psu.edu/cab38/ColorBrewer/ColorBrewer_updates.html +% +% +% [colormap]=cbrewer(ctype, cname, ncol, interp_method) +% +% INPUT: +% - ctype: type of color table 'seq' (sequential), 'div' (diverging), 'qual' (qualitative) +% - cname: name of colortable. It changes depending on ctype. +% - ncol: number of color in the table. It changes according to ctype and +% cname +% - interp_method: interpolation method (see interp1.m). Default is "cubic" ) +% +% A note on the number of colors: Based on the original data, there is +% only a certain number of colors available for each type and name of +% colortable. When 'ncol' is larger then the maximum number of colors +% originally given, an interpolation routine is called (interp1) to produce +% the "extended" colormaps. +% +% Example: To produce a colortable CT of ncol X 3 entries (RGB) of +% sequential type and named 'Blues' with 8 colors: +% CT=cbrewer('seq', 'Blues', 8); +% To use this colortable as colormap, simply call: +% colormap(CT) +% +% To see the various colormaps available according to their types and +% names, simply call: cbrewer() +% +% This product includes color specifications and designs developed by +% Cynthia Brewer (http://colorbrewer.org/). +% +% Author: Charles Robert +% email: tannoudji@hotmail.com +% Date: 06.12.2011 +% ------------------------------ +% 18.09.2015 Minor fixes, fixed a bug where the 'spectral' color table did not appear in the preview + + +% load colorbrewer data +load('colorbrewer.mat') +% initialise the colormap is there are any problems +colormap=[]; +if (~exist('interp_method', 'var')) + interp_method='cubic'; +end + +% If no arguments +if (~exist('ctype', 'var') | ~exist('cname', 'var') | ~exist('ncol', 'var')) + disp(' ') + disp('[colormap] = cbrewer(ctype, cname, ncol [, interp_method])') + disp(' ') + disp('INPUT:') + disp(' - ctype: type of color table *seq* (sequential), *div* (divergent), *qual* (qualitative)') + disp(' - cname: name of colortable. It changes depending on ctype.') + disp(' - ncol: number of color in the table. It changes according to ctype and cname') + disp(' - interp_method: interpolation method (see interp1.m). Default is "cubic" )') + + disp(' ') + disp('Sequential tables:') + z={'Blues','BuGn','BuPu','GnBu','Greens','Greys','Oranges','OrRd','PuBu','PuBuGn','PuRd',... + 'Purples','RdPu', 'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd', 'Spectral'}; + disp(z') + + disp('Divergent tables:') + z={'BrBG', 'PiYG', 'PRGn', 'PuOr', 'RdBu', 'RdGy', 'RdYlBu', 'RdYlGn'}; + disp(z') + + disp(' ') + disp('Qualitative tables:') + %getfield(colorbrewer, 'qual') + z={'Accent', 'Dark2', 'Paired', 'Pastel1', 'Pastel2', 'Set1', 'Set2', 'Set3'}; + disp(z') + + plot_brewer_cmap + return +end + +% Verify that the input is appropriate +ctype_names={'div', 'seq', 'qual'}; +if (~ismember(ctype,ctype_names)) + disp('ctype must be either: *div*, *seq* or *qual*') + colormap=[]; + return +end + +if (~isfield(colorbrewer.(ctype),cname)) + disp(['The name of the colortable of type *' ctype '* must be one of the following:']) + getfield(colorbrewer, ctype) + colormap=[]; + return +end + +if (ncol>length(colorbrewer.(ctype).(cname))) +% disp(' ') +% disp('----------------------------------------------------------------------') +% disp(['The maximum number of colors for table *' cname '* is ' num2str(length(colorbrewer.(ctype).(cname)))]) +% disp(['The new colormap will be extrapolated from these ' num2str(length(colorbrewer.(ctype).(cname))) ' values']) +% disp('----------------------------------------------------------------------') +% disp(' ') + cbrew_init=colorbrewer.(ctype).(cname){length(colorbrewer.(ctype).(cname))}; + colormap=interpolate_cbrewer(cbrew_init, interp_method, ncol); + colormap=colormap./255; + return +end + +if (isempty(colorbrewer.(ctype).(cname){ncol})) + + while(isempty(colorbrewer.(ctype).(cname){ncol})) + ncol=ncol+1; + end + disp(' ') + disp('----------------------------------------------------------------------') + disp(['The minimum number of colors for table *' cname '* is ' num2str(ncol)]) + disp('This minimum value shall be defined as ncol instead') + disp('----------------------------------------------------------------------') + disp(' ') +end + +colormap=(colorbrewer.(ctype).(cname){ncol})./255; + +end \ No newline at end of file diff --git a/dependencies/cbrewer/cbrewer/cbrewer_preview.jpg b/dependencies/cbrewer/cbrewer/cbrewer_preview.jpg new file mode 100644 index 0000000..bd2830a Binary files /dev/null and b/dependencies/cbrewer/cbrewer/cbrewer_preview.jpg differ diff --git a/dependencies/cbrewer/cbrewer/change_jet.m b/dependencies/cbrewer/cbrewer/change_jet.m new file mode 100644 index 0000000..b8d4ecb --- /dev/null +++ b/dependencies/cbrewer/cbrewer/change_jet.m @@ -0,0 +1,64 @@ +% This script help produce a new 'jet'-like colormap based on other RGB reference colors + +% ------- I WAS ASKED --------------- +% "is there a chance that you could add a diverging map going from blue to green to red as in jet, +% but using the red and blue from your RdBu map and the third darkest green from your RdYlGn map?"" +% +% ANSWER: +% You should construct the new colormap based on the existing RGB values of 'jet' +% but projecting these RGB values on your new RGB basis. +% ----------------------------------- + +% load colormaps +jet=colormap('jet'); +RdBu=cbrewer('div', 'RdBu', 11); +RdYlGn=cbrewer('div', 'RdYlGn', 11); + +% Define the new R, G, B references (p stands for prime) +Rp=RdBu(1,:); +Bp=RdBu(end, :); +Gp=RdYlGn(end-2, :); +RGBp=[Rp;Gp;Bp]; + +% construct the new colormap based on the existing RGB values of jet +% Project the RGB values on your new basis +newjet = jet*RGBp; + +% store data in a strcuture, easier to handle +cmap.jet=jet; +cmap.newjet=newjet; +cnames={'jet', 'newjet'}; + +% plot the RGB values +fh=figure(); +colors={'r', 'g', 'b'}; +for iname=1:length(cnames) + subplot(length(cnames),1,iname) + dat=cmap.(cnames{end-iname+1}); + for icol=1:size(dat,2) + plot(dat(:,icol), 'color', colors{icol}, 'linewidth', 2);hold on; + end % icol + title([' "' cnames{end-iname+1} '" in RGB plot']) +end + +% plot the colormaps +fh=figure(); +for iname=1:length(cnames) + F=cmap.(cnames{iname}); + ncol=length(F); + fg=1./ncol; % geometrical factor + X=fg.*[0 0 1 1]; + Y=0.1.*[1 0 0 1]+(2*iname-1)*0.1; + + for icol=1:ncol + X2=X+fg.*(icol-1); + fill(X2,Y,F(icol, :), 'linestyle', 'none') + hold all + end % icol + text(-0.1, mean(Y), cnames{iname}, 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'FontSize', 10, 'FontName' , 'AvantGarde') + xlim([-0.4, 1]) + axis off + set(gcf, 'color', [1 1 1]) + ylim([0.1 1.05.*max(Y)]); + end % iname + diff --git a/dependencies/cbrewer/cbrewer/colorbrewer.mat b/dependencies/cbrewer/cbrewer/colorbrewer.mat new file mode 100644 index 0000000..ec59ef4 Binary files /dev/null and b/dependencies/cbrewer/cbrewer/colorbrewer.mat differ diff --git a/dependencies/cbrewer/cbrewer/interpolate_cbrewer.m b/dependencies/cbrewer/cbrewer/interpolate_cbrewer.m new file mode 100644 index 0000000..e8b5e21 --- /dev/null +++ b/dependencies/cbrewer/cbrewer/interpolate_cbrewer.m @@ -0,0 +1,36 @@ +function [interp_cmap]=interpolate_cbrewer(cbrew_init, interp_method, ncolors) +% +% INTERPOLATE_CBREWER - interpolate a colorbrewer map to ncolors levels +% +% INPUT: +% - cbrew_init: the initial colormap with format N*3 +% - interp_method: interpolation method, which can be the following: +% 'nearest' - nearest neighbor interpolation +% 'linear' - bilinear interpolation +% 'spline' - spline interpolation +% 'cubic' - bicubic interpolation as long as the data is +% uniformly spaced, otherwise the same as 'spline' +% - ncolors=desired number of colors +% +% Author: Charles Robert +% email: tannoudji@hotmail.com +% Date: 14.10.2011 + + +% just to make sure, in case someone puts in a decimal +ncolors=round(ncolors); + +% How many data points of the colormap available +nmax=size(cbrew_init,1); + +% create the associated X axis (using round to get rid of decimals) +a=(ncolors-1)./(nmax-1); +X=round([0 a:a:(ncolors-1)]); +X2=0:ncolors-1; + +z=interp1(X,cbrew_init(:,1),X2,interp_method); +z2=interp1(X,cbrew_init(:,2),X2,interp_method); +z3=interp1(X,cbrew_init(:,3),X2, interp_method); +interp_cmap=round([z' z2' z3']); + +end \ No newline at end of file diff --git a/dependencies/cbrewer/cbrewer/plot_brewer_cmap.m b/dependencies/cbrewer/cbrewer/plot_brewer_cmap.m new file mode 100644 index 0000000..a5cab9e --- /dev/null +++ b/dependencies/cbrewer/cbrewer/plot_brewer_cmap.m @@ -0,0 +1,50 @@ +% Plots and identifies the various colorbrewer tables available. +% Is called by cbrewer.m when no arguments are given. +% +% Author: Charles Robert +% email: tannoudji@hotmail.com +% Date: 14.10.2011 + + + +load('colorbrewer.mat') + +ctypes={'div', 'seq', 'qual'}; +ctypes_title={'Diverging', 'Sequential', 'Qualitative'}; +cnames{1,:}={'BrBG', 'PiYG', 'PRGn', 'PuOr', 'RdBu', 'RdGy', 'RdYlBu', 'RdYlGn', 'Spectral'}; +cnames{2,:}={'Blues','BuGn','BuPu','GnBu','Greens','Greys','Oranges','OrRd','PuBu','PuBuGn','PuRd',... + 'Purples','RdPu', 'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd'}; +cnames{3,:}={'Accent', 'Dark2', 'Paired', 'Pastel1', 'Pastel2', 'Set1', 'Set2', 'Set3'}; + +figure('position', [314 327 807 420]) +for itype=1:3 + + %fh(itype)=figure(); + subplot(1,3,itype) + + for iname=1:length(cnames{itype,:}) + + ncol=length(colorbrewer.(ctypes{itype}).(cnames{itype}{iname})); + fg=1./ncol; % geometrical factor + + X=fg.*[0 0 1 1]; + Y=0.1.*[1 0 0 1]+(2*iname-1)*0.1; + F=cbrewer(ctypes{itype}, cnames{itype}{iname}, ncol); + + for icol=1:ncol + X2=X+fg.*(icol-1); + fill(X2,Y,F(icol, :), 'linestyle', 'none') + text(-0.1, mean(Y), cnames{itype}{iname}, 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'FontSize', 10, 'FontName' , 'AvantGarde') + xlim([-0.4, 1]) + hold all + end % icol + %set(gca, 'box', 'off') + title(ctypes_title{itype}, 'FontWeight', 'bold', 'FontSize', 16, 'FontName' , 'AvantGarde') + axis off + set(gcf, 'color', [1 1 1]) + end % iname + ylim([0.1 1.05.*max(Y)]); +end %itype + +set(gcf, 'MenuBar', 'none') +set(gcf, 'Name', 'ColorBrewer Color maps') \ No newline at end of file diff --git a/dependencies/cbrewer/license.txt b/dependencies/cbrewer/license.txt new file mode 100644 index 0000000..ded8179 --- /dev/null +++ b/dependencies/cbrewer/license.txt @@ -0,0 +1,24 @@ +Copyright (c) 2011, Charles Robert +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the distribution + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. diff --git a/dependencies/make_nii.m b/dependencies/make_nii.m new file mode 100644 index 0000000..1af1c5e --- /dev/null +++ b/dependencies/make_nii.m @@ -0,0 +1,256 @@ +% Make NIfTI structure specified by an N-D matrix. Usually, N is 3 for +% 3D matrix [x y z], or 4 for 4D matrix with time series [x y z t]. +% Optional parameters can also be included, such as: voxel_size, +% origin, datatype, and description. +% +% Once the NIfTI structure is made, it can be saved into NIfTI file +% using "save_nii" command (for more detail, type: help save_nii). +% +% Usage: nii = make_nii(img, [voxel_size], [origin], [datatype], [description]) +% +% Where: +% +% img: Usually, img is a 3D matrix [x y z], or a 4D +% matrix with time series [x y z t]. However, +% NIfTI allows a maximum of 7D matrix. When the +% image is in RGB format, make sure that the size +% of 4th dimension is always 3 (i.e. [R G B]). In +% that case, make sure that you must specify RGB +% datatype, which is either 128 or 511. +% +% voxel_size (optional): Voxel size in millimeter for each +% dimension. Default is [1 1 1]. +% +% origin (optional): The AC origin. Default is [0 0 0]. +% +% datatype (optional): Storage data type: +% 2 - uint8, 4 - int16, 8 - int32, 16 - float32, +% 32 - complex64, 64 - float64, 128 - RGB24, +% 256 - int8, 511 - RGB96, 512 - uint16, +% 768 - uint32, 1792 - complex128 +% Default will use the data type of 'img' matrix +% For RGB image, you must specify it to either 128 +% or 511. +% +% description (optional): Description of data. Default is ''. +% +% e.g.: +% origin = [33 44 13]; datatype = 64; +% nii = make_nii(img, [], origin, datatype); % default voxel_size +% +% NIFTI data format can be found on: http://nifti.nimh.nih.gov +% +% - Jimmy Shen (jimmy@rotman-baycrest.on.ca) +% +function nii = make_nii(varargin) + + nii.img = varargin{1}; + dims = size(nii.img); + dims = [length(dims) dims ones(1,8)]; + dims = dims(1:8); + + voxel_size = [0 ones(1,7)]; + origin = zeros(1,5); + descrip = ''; + + switch class(nii.img) + case 'uint8' + datatype = 2; + case 'int16' + datatype = 4; + case 'int32' + datatype = 8; + case 'single' + if isreal(nii.img) + datatype = 16; + else + datatype = 32; + end + case 'double' + if isreal(nii.img) + datatype = 64; + else + datatype = 1792; + end + case 'int8' + datatype = 256; + case 'uint16' + datatype = 512; + case 'uint32' + datatype = 768; + otherwise + error('Datatype is not supported by make_nii.'); + end + + if nargin > 1 & ~isempty(varargin{2}) + voxel_size(2:4) = double(varargin{2}); + end + + if nargin > 2 & ~isempty(varargin{3}) + origin(1:3) = double(varargin{3}); + end + + if nargin > 3 & ~isempty(varargin{4}) + datatype = double(varargin{4}); + + if datatype == 128 | datatype == 511 + dims(5) = []; + dims(1) = dims(1) - 1; + dims = [dims 1]; + end + end + + if nargin > 4 & ~isempty(varargin{5}) + descrip = varargin{5}; + end + + if ndims(nii.img) > 7 + error('NIfTI only allows a maximum of 7 Dimension matrix.'); + end + + maxval = round(double(max(nii.img(:)))); + minval = round(double(min(nii.img(:)))); + + nii.hdr = make_header(dims, voxel_size, origin, datatype, ... + descrip, maxval, minval); + + switch nii.hdr.dime.datatype + case 2 + nii.img = uint8(nii.img); + case 4 + nii.img = int16(nii.img); + case 8 + nii.img = int32(nii.img); + case 16 + nii.img = single(nii.img); + case 32 + nii.img = single(nii.img); + case 64 + nii.img = double(nii.img); + case 128 + nii.img = uint8(nii.img); + case 256 + nii.img = int8(nii.img); + case 511 + img = double(nii.img(:)); + img = single((img - min(img))/(max(img) - min(img))); + nii.img = reshape(img, size(nii.img)); + nii.hdr.dime.glmax = double(max(img)); + nii.hdr.dime.glmin = double(min(img)); + case 512 + nii.img = uint16(nii.img); + case 768 + nii.img = uint32(nii.img); + case 1792 + nii.img = double(nii.img); + otherwise + error('Datatype is not supported by make_nii.'); + end + + return; % make_nii + + +%--------------------------------------------------------------------- +function hdr = make_header(dims, voxel_size, origin, datatype, ... + descrip, maxval, minval) + + hdr.hk = header_key; + hdr.dime = image_dimension(dims, voxel_size, datatype, maxval, minval); + hdr.hist = data_history(origin, descrip); + + return; % make_header + + +%--------------------------------------------------------------------- +function hk = header_key + + hk.sizeof_hdr = 348; % must be 348! + hk.data_type = ''; + hk.db_name = ''; + hk.extents = 0; + hk.session_error = 0; + hk.regular = 'r'; + hk.dim_info = 0; + + return; % header_key + + +%--------------------------------------------------------------------- +function dime = image_dimension(dims, voxel_size, datatype, maxval, minval) + + dime.dim = dims; + dime.intent_p1 = 0; + dime.intent_p2 = 0; + dime.intent_p3 = 0; + dime.intent_code = 0; + dime.datatype = datatype; + + switch dime.datatype + case 2, + dime.bitpix = 8; precision = 'uint8'; + case 4, + dime.bitpix = 16; precision = 'int16'; + case 8, + dime.bitpix = 32; precision = 'int32'; + case 16, + dime.bitpix = 32; precision = 'float32'; + case 32, + dime.bitpix = 64; precision = 'float32'; + case 64, + dime.bitpix = 64; precision = 'float64'; + case 128 + dime.bitpix = 24; precision = 'uint8'; + case 256 + dime.bitpix = 8; precision = 'int8'; + case 511 + dime.bitpix = 96; precision = 'float32'; + case 512 + dime.bitpix = 16; precision = 'uint16'; + case 768 + dime.bitpix = 32; precision = 'uint32'; + case 1792, + dime.bitpix = 128; precision = 'float64'; + otherwise + error('Datatype is not supported by make_nii.'); + end + + dime.slice_start = 0; + dime.pixdim = voxel_size; + dime.vox_offset = 0; + dime.scl_slope = 0; + dime.scl_inter = 0; + dime.slice_end = 0; + dime.slice_code = 0; + dime.xyzt_units = 0; + dime.cal_max = 0; + dime.cal_min = 0; + dime.slice_duration = 0; + dime.toffset = 0; + dime.glmax = maxval; + dime.glmin = minval; + + return; % image_dimension + + +%--------------------------------------------------------------------- +function hist = data_history(origin, descrip) + + hist.descrip = descrip; + hist.aux_file = 'none'; + hist.qform_code = 0; + hist.sform_code = 0; + hist.quatern_b = 0; + hist.quatern_c = 0; + hist.quatern_d = 0; + hist.qoffset_x = 0; + hist.qoffset_y = 0; + hist.qoffset_z = 0; + hist.srow_x = zeros(1,4); + hist.srow_y = zeros(1,4); + hist.srow_z = zeros(1,4); + hist.intent_name = ''; + hist.magic = ''; + hist.originator = origin; + + return; % data_history + diff --git a/dependencies/save_nii.m b/dependencies/save_nii.m new file mode 100644 index 0000000..b2d0dd4 --- /dev/null +++ b/dependencies/save_nii.m @@ -0,0 +1,286 @@ +% Save NIFTI dataset. Support both *.nii and *.hdr/*.img file extension. +% If file extension is not provided, *.hdr/*.img will be used as default. +% +% Usage: save_nii(nii, filename, [old_RGB]) +% +% nii.hdr - struct with NIFTI header fields (from load_nii.m or make_nii.m) +% +% nii.img - 3D (or 4D) matrix of NIFTI data. +% +% filename - NIFTI file name. +% +% old_RGB - an optional boolean variable to handle special RGB data +% sequence [R1 R2 ... G1 G2 ... B1 B2 ...] that is used only by +% AnalyzeDirect (Analyze Software). Since both NIfTI and Analyze +% file format use RGB triple [R1 G1 B1 R2 G2 B2 ...] sequentially +% for each voxel, this variable is set to FALSE by default. If you +% would like the saved image only to be opened by AnalyzeDirect +% Software, set old_RGB to TRUE (or 1). It will be set to 0, if it +% is default or empty. +% +% Tip: to change the data type, set nii.hdr.dime.datatype, +% and nii.hdr.dime.bitpix to: +% +% 0 None (Unknown bit per voxel) % DT_NONE, DT_UNKNOWN +% 1 Binary (ubit1, bitpix=1) % DT_BINARY +% 2 Unsigned char (uchar or uint8, bitpix=8) % DT_UINT8, NIFTI_TYPE_UINT8 +% 4 Signed short (int16, bitpix=16) % DT_INT16, NIFTI_TYPE_INT16 +% 8 Signed integer (int32, bitpix=32) % DT_INT32, NIFTI_TYPE_INT32 +% 16 Floating point (single or float32, bitpix=32) % DT_FLOAT32, NIFTI_TYPE_FLOAT32 +% 32 Complex, 2 float32 (Use float32, bitpix=64) % DT_COMPLEX64, NIFTI_TYPE_COMPLEX64 +% 64 Double precision (double or float64, bitpix=64) % DT_FLOAT64, NIFTI_TYPE_FLOAT64 +% 128 uint RGB (Use uint8, bitpix=24) % DT_RGB24, NIFTI_TYPE_RGB24 +% 256 Signed char (schar or int8, bitpix=8) % DT_INT8, NIFTI_TYPE_INT8 +% 511 Single RGB (Use float32, bitpix=96) % DT_RGB96, NIFTI_TYPE_RGB96 +% 512 Unsigned short (uint16, bitpix=16) % DT_UNINT16, NIFTI_TYPE_UNINT16 +% 768 Unsigned integer (uint32, bitpix=32) % DT_UNINT32, NIFTI_TYPE_UNINT32 +% 1024 Signed long long (int64, bitpix=64) % DT_INT64, NIFTI_TYPE_INT64 +% 1280 Unsigned long long (uint64, bitpix=64) % DT_UINT64, NIFTI_TYPE_UINT64 +% 1536 Long double, float128 (Unsupported, bitpix=128) % DT_FLOAT128, NIFTI_TYPE_FLOAT128 +% 1792 Complex128, 2 float64 (Use float64, bitpix=128) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX128 +% 2048 Complex256, 2 float128 (Unsupported, bitpix=256) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX128 +% +% Part of this file is copied and modified from: +% http://www.mathworks.com/matlabcentral/fileexchange/1878-mri-analyze-tools +% +% NIFTI data format can be found on: http://nifti.nimh.nih.gov +% +% - Jimmy Shen (jimmy@rotman-baycrest.on.ca) +% - "old_RGB" related codes in "save_nii.m" are added by Mike Harms (2006.06.28) +% +function save_nii(nii, fileprefix, old_RGB) + + if ~exist('nii','var') | isempty(nii) | ~isfield(nii,'hdr') | ... + ~isfield(nii,'img') | ~exist('fileprefix','var') | isempty(fileprefix) + + error('Usage: save_nii(nii, filename, [old_RGB])'); + end + + if isfield(nii,'untouch') & nii.untouch == 1 + error('Usage: please use ''save_untouch_nii.m'' for the untouched structure.'); + end + + if ~exist('old_RGB','var') | isempty(old_RGB) + old_RGB = 0; + end + + v = version; + + % Check file extension. If .gz, unpack it into temp folder + % + if length(fileprefix) > 2 & strcmp(fileprefix(end-2:end), '.gz') + + if ~strcmp(fileprefix(end-6:end), '.img.gz') & ... + ~strcmp(fileprefix(end-6:end), '.hdr.gz') & ... + ~strcmp(fileprefix(end-6:end), '.nii.gz') + + error('Please check filename.'); + end + + if str2num(v(1:3)) < 7.1 | ~usejava('jvm') + error('Please use MATLAB 7.1 (with java) and above, or run gunzip outside MATLAB.'); + else + gzFile = 1; + fileprefix = fileprefix(1:end-3); + end + end + + filetype = 1; + + % Note: fileprefix is actually the filename you want to save + % + if findstr('.nii',fileprefix) & strcmp(fileprefix(end-3:end), '.nii') + filetype = 2; + fileprefix(end-3:end)=''; + end + + if findstr('.hdr',fileprefix) & strcmp(fileprefix(end-3:end), '.hdr') + fileprefix(end-3:end)=''; + end + + if findstr('.img',fileprefix) & strcmp(fileprefix(end-3:end), '.img') + fileprefix(end-3:end)=''; + end + + write_nii(nii, filetype, fileprefix, old_RGB); + + % gzip output file if requested + % + if exist('gzFile', 'var') + if filetype == 1 + gzip([fileprefix, '.img']); + delete([fileprefix, '.img']); + gzip([fileprefix, '.hdr']); + delete([fileprefix, '.hdr']); + elseif filetype == 2 + gzip([fileprefix, '.nii']); + delete([fileprefix, '.nii']); + end; + end; + + if filetype == 1 + + % So earlier versions of SPM can also open it with correct originator + % + M=[[diag(nii.hdr.dime.pixdim(2:4)) -[nii.hdr.hist.originator(1:3).*nii.hdr.dime.pixdim(2:4)]'];[0 0 0 1]]; + save([fileprefix '.mat'], 'M'); + end + + return % save_nii + + +%----------------------------------------------------------------------------------- +function write_nii(nii, filetype, fileprefix, old_RGB) + + hdr = nii.hdr; + + if isfield(nii,'ext') & ~isempty(nii.ext) + ext = nii.ext; + [ext, esize_total] = verify_nii_ext(ext); + else + ext = []; + end + + switch double(hdr.dime.datatype), + case 1, + hdr.dime.bitpix = int16(1 ); precision = 'ubit1'; + case 2, + hdr.dime.bitpix = int16(8 ); precision = 'uint8'; + case 4, + hdr.dime.bitpix = int16(16); precision = 'int16'; + case 8, + hdr.dime.bitpix = int16(32); precision = 'int32'; + case 16, + hdr.dime.bitpix = int16(32); precision = 'float32'; + case 32, + hdr.dime.bitpix = int16(64); precision = 'float32'; + case 64, + hdr.dime.bitpix = int16(64); precision = 'float64'; + case 128, + hdr.dime.bitpix = int16(24); precision = 'uint8'; + case 256 + hdr.dime.bitpix = int16(8 ); precision = 'int8'; + case 511, + hdr.dime.bitpix = int16(96); precision = 'float32'; + case 512 + hdr.dime.bitpix = int16(16); precision = 'uint16'; + case 768 + hdr.dime.bitpix = int16(32); precision = 'uint32'; + case 1024 + hdr.dime.bitpix = int16(64); precision = 'int64'; + case 1280 + hdr.dime.bitpix = int16(64); precision = 'uint64'; + case 1792, + hdr.dime.bitpix = int16(128); precision = 'float64'; + otherwise + error('This datatype is not supported'); + end + + hdr.dime.glmax = round(double(max(nii.img(:)))); + hdr.dime.glmin = round(double(min(nii.img(:)))); + + if filetype == 2 + fid = fopen(sprintf('%s.nii',fileprefix),'w'); + + if fid < 0, + msg = sprintf('Cannot open file %s.nii.',fileprefix); + error(msg); + end + + hdr.dime.vox_offset = 352; + + if ~isempty(ext) + hdr.dime.vox_offset = hdr.dime.vox_offset + esize_total; + end + + hdr.hist.magic = 'n+1'; + save_nii_hdr(hdr, fid); + + if ~isempty(ext) + save_nii_ext(ext, fid); + end + else + fid = fopen(sprintf('%s.hdr',fileprefix),'w'); + + if fid < 0, + msg = sprintf('Cannot open file %s.hdr.',fileprefix); + error(msg); + end + + hdr.dime.vox_offset = 0; + hdr.hist.magic = 'ni1'; + save_nii_hdr(hdr, fid); + + if ~isempty(ext) + save_nii_ext(ext, fid); + end + + fclose(fid); + fid = fopen(sprintf('%s.img',fileprefix),'w'); + end + + ScanDim = double(hdr.dime.dim(5)); % t + SliceDim = double(hdr.dime.dim(4)); % z + RowDim = double(hdr.dime.dim(3)); % y + PixelDim = double(hdr.dime.dim(2)); % x + SliceSz = double(hdr.dime.pixdim(4)); + RowSz = double(hdr.dime.pixdim(3)); + PixelSz = double(hdr.dime.pixdim(2)); + + x = 1:PixelDim; + + if filetype == 2 & isempty(ext) + skip_bytes = double(hdr.dime.vox_offset) - 348; + else + skip_bytes = 0; + end + + if double(hdr.dime.datatype) == 128 + + % RGB planes are expected to be in the 4th dimension of nii.img + % + if(size(nii.img,4)~=3) + error(['The NII structure does not appear to have 3 RGB color planes in the 4th dimension']); + end + + if old_RGB + nii.img = permute(nii.img, [1 2 4 3 5 6 7 8]); + else + nii.img = permute(nii.img, [4 1 2 3 5 6 7 8]); + end + end + + if double(hdr.dime.datatype) == 511 + + % RGB planes are expected to be in the 4th dimension of nii.img + % + if(size(nii.img,4)~=3) + error(['The NII structure does not appear to have 3 RGB color planes in the 4th dimension']); + end + + if old_RGB + nii.img = permute(nii.img, [1 2 4 3 5 6 7 8]); + else + nii.img = permute(nii.img, [4 1 2 3 5 6 7 8]); + end + end + + % For complex float32 or complex float64, voxel values + % include [real, imag] + % + if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792 + real_img = real(nii.img(:))'; + nii.img = imag(nii.img(:))'; + nii.img = [real_img; nii.img]; + end + + if skip_bytes + fwrite(fid, zeros(1,skip_bytes), 'uint8'); + end + + fwrite(fid, nii.img, precision); +% fwrite(fid, nii.img, precision, skip_bytes); % error using skip + fclose(fid); + + return; % write_nii + diff --git a/ppicaps/deconvolveSubj.m b/ppicaps/deconvolveSubj.m new file mode 100644 index 0000000..9d7550d --- /dev/null +++ b/ppicaps/deconvolveSubj.m @@ -0,0 +1,130 @@ +% Function copied from the relevant section of SPM8's spm_peb_ppi.m +% function for devonvolution, usually used in the context of PPI analyses. +% +% Last checked: September 2019 + +function deconvData = deconvolveSubj(subj) + + +% Load SPM +%-------------------------------------------------------------------------- +analyze_dir = [subj.dataDir 'GLM_SciFun/']; +SPM = spm_select('FPListRec',analyze_dir,'SPM.mat$'); load(SPM); +tempVOI = load([subj.dataDir '/VOI_PCC_1.mat']); +xY = tempVOI.xY; +VI = SPM.xY.VY; + +% Setup variables +%-------------------------------------------------------------------------- +RT = SPM.xY.RT; +dt = SPM.xBF.dt; +NT = round(RT/dt); +N = length(xY(1).u); +k = 1:NT:N*NT; +Sess = SPM.Sess(1); + +% Create basis functions and hrf in scan time and microtime +%-------------------------------------------------------------------------- +hrf = spm_hrf(dt); % 1 corresponds to the microtime resolution I want + + +% Create convolved explanatory {Hxb} variables in scan time +%-------------------------------------------------------------------------- +xb = spm_dctmtx(N*NT + 128,N); %xb = DCT matrix +Hxb = zeros(N,N); +for i = 1:N + Hx = conv(xb(:,i),hrf); + Hxb(:,i) = Hx(k + 128); +end +xb = xb(129:end,:); + +% Get confounds (in scan time) and constant term +%-------------------------------------------------------------------------- +X0 = xY(1).X0; +M = size(X0,2); + +% Specify covariance components; assume neuronal response is white +% treating confounds as fixed effects +%-------------------------------------------------------------------------- +Q = speye(N,N)*N/trace(Hxb'*Hxb); +Q = blkdiag(Q, speye(M,M)*1e6 ); + + +% Get whitening matrix (NB: confounds have already been whitened) +%-------------------------------------------------------------------------- +W = SPM.xX.W(Sess.row,Sess.row); + + +% Create structure for spm_PEB +%-------------------------------------------------------------------------- +clear P +P{1}.X = [W*Hxb X0]; % Design matrix for lowest level +P{1}.C = speye(N,N)/4; % i.i.d assumptions +P{2}.X = sparse(N + M,1); % Design matrix for parameters (0's) +P{2}.C = Q; + +% Deconvolution +%========================================================================= +tic +% Iterate over voxels and deconvolve +total = 153594; counter = 0; percentage = 0; +x= SPM.xY.VY(1).dim(1); +y= SPM.xY.VY(1).dim(2); +z= SPM.xY.VY(1).dim(3); +t = SPM.nscan; +deconvData = zeros(x,y,z,t); + +% cd into analyze_dir to be able to calculate beta +cd(analyze_dir); +for ix = 1:x + for iy = 1:y + for iz = 1:z + + thisVoxel = spm_data_read(VI,'xyz',[ix;iy;iz]); + thisVoxel = spm_filter(SPM.xX.K,W*thisVoxel); + + % Remove null space of contrast + beta = spm_data_read(SPM.Vbeta,'xyz',[ix;iy;iz]); + if ~isnan(mean(beta)) + thisVoxel = thisVoxel - spm_FcUtil('Y0',SPM.xCon(xY.Ic),SPM.xX.xKXs,beta); + end + % Simple deconvolution + % ------------------------------------------------------------- + C = spm_PEB(thisVoxel,P); + xn = xb*C{2}.E(1:N); + xn = spm_detrend(xn); + + % Save variables (NOTE: xn is in microtime and does not account for + % slice timing shifts). To convert to BOLD signal convolve with a hrf. + % Use a microtime to scan time index to convert to scan time: e.g., + % k = 1:NT:N*NT; where NT = number of bins per TR = TR/dt or SPM.xBF.T + % and N = number of scans in the session. Finally account for slice + % timing effects by shifting the index accordingly. + %---------------------------------------------------------------------- + deconvVoxel = downsample(xn,16); + deconvData(ix, iy, iz,:) = deconvVoxel; + + counter = counter +1; + percentage = (counter / total ) * 100; + + end + end +end +toc + +% CD back into working folder to save file +cd(pwd); +save([subj.curSubj '_deconvolvedData.mat'], 'deconvData'); + + + +end + + + + + + + + + diff --git a/ppicaps/initialize_vars.m b/ppicaps/initialize_vars.m new file mode 100644 index 0000000..f3aa8bd --- /dev/null +++ b/ppicaps/initialize_vars.m @@ -0,0 +1,37 @@ +% This function initialises the subject variables for the batch analysis of +% the Movie Watching data. +% _________________________________________________________________________ +% Last checked: September 2019 + +% Lorena Freitas +% $Id: initialize_vars.m 11 2016-12-01 16:26:24F Lorena $ + + +function [b] = initialize_vars(subject,group) + +% SPM info +%-------------------------------------------------------------------------- +b.spmDir = fileparts(which('spm')); % path to SPM installation + + +% Directory information +%-------------------------------------------------------------------------- +dataDir = ''; + + +% Subject information +%-------------------------------------------------------------------------- +b.group = group; +b.curSubj = subject; % subject code (e.g., 'sujet01') +b.dataDir = strcat(dataDir, filesep , b.curSubj, filesep); % make data directory subject-specific + +% Data type to use (deconvolved = 1 on non-deconvolved = 0) +%-------------------------------------------------------------------------- +b.deconvolveData = 1; + +% Folder for functional files +%-------------------------------------------------------------------------- +b.dir_fonc = ''; + + +end diff --git a/ppicaps/kmeanspp.m b/ppicaps/kmeanspp.m new file mode 100644 index 0000000..0bb5623 --- /dev/null +++ b/ppicaps/kmeanspp.m @@ -0,0 +1,179 @@ +function [L,Ix,C,dis,sD] = kmeanspp(X,k, dist,varargin) +% KMEANSPP ? Cluster multivariate data using the k-means++ algorithm. +% [L,Ix,C,dst,sD] = kmeans(X,k) takes an n-by-p matrix X as input containing p observations of +% n-dimensional data and produces a 1-by-size(X,2) vector L with one class +% label per column in X, a size(X,1)-by-k matrix C containing the +% centers corresponding to each class, a vector Ix indicating the +% sign of the time point (1 = positive; 2 = negative) a vector dis +% containing the average distance from each centroid to their corresponding +% points and a vector sD indictating the standard deviation of each +% cluster +% +% Author: Lorena Freitas (Lorena.Freitas@epfl.ch) +% +% % Last checked: September 2019 +% +% 'varargin' +% #replicates +% or C - matrix of centroids +% +% References: +% [1] J. B. MacQueen, "Some Methods for Classification and Analysis of +% MultiVariate Observations", in Proc. of the fifth Berkeley +% Symposium on Mathematical Statistics and Probability, L. M. L. Cam +% and J. Neyman, eds., vol. 1, UC Press, 1967, pp. 281-297. +% [2] D. Arthur and S. Vassilvitskii, "k-means++: The Advantages of +% Careful Seeding", Technical Report 2006-13, Stanford InfoLab, 2006. + + +% reseed random number generator +%rng shuffle + +% Default value of replicates: only used when centroids are passed as an +% argument, which means we only need to assign each frame to a centroid, +% and therefore there's no need to run it more than once. +replicates = 1; + +switch dist + case 'sqEuclidean' + scoremethod = @sqEuclideanScore; + distancemethod = @sqEuclideanDistance; + case 'cosine' + scoremethod = @CosineScore; + distancemethod = @CosineDistance; + + % Normalise vectors to unit norm + X = X./repmat(sqrt(sum(X.*X)),size(X,1),1); +end + +if nargin > 3 + + % If 3rd argument is a matrix of centroids, assign frames to centroids + if length(varargin{1}) > 1 + C = varargin{1}; + [M1,Lpos] = min(scoremethod(C,X),[],1); + [M2,Lneg] = min(scoremethod(-C,X),[],1); + [~,Ix] = min([M1;M2]); + L = ones(1,length(Lpos)); L(Ix==1) = Lpos(Ix==1); L(Ix==2) = Lneg(Ix==2); + M = ones(1,length(M1)); M(Ix==1) = M1(Ix==1); M(Ix==2) = M2(Ix==2); + + % if 3rd argument is the number of Ks + else + replicates = varargin{1}; + bestSolution.optVar = Inf; %realmin; + + for repl=1:replicates + + L = []; + L1 = 0; + + while length(unique(L)) ~= k + + % The k-means++ initialization. The code below is inspired + % by Laurent Sorber's code at MATLAB Central File Exchange + % from 2013-02-08. + C = X(:,1+round(rand*(size(X,2)-1))); + L = ones(1,size(X,2)); + for i = 2:k + D = X-C(:,L); D2 = X+C(:,L); % D2 is distance to negative centroid; + D = min([cumsum(sqrt(dot(D,D,1)));cumsum(sqrt(dot(D2,D2,1)))]); + if D(end) == 0, C(:,i:k) = X(:,ones(1,k-i+1)); return; end + C(:,i) = X(:,find(rand < D/D(end),1)); + + [M1,Lpos] = min(scoremethod(X,C)); + [M2,Lneg] = min(scoremethod(X,-C)); + [~,Ix] = min([M1;M2]); + L = ones(1,length(Lpos)); L(Ix==1) = Lpos(Ix==1); L(Ix==2) = Lneg(Ix==2); + end + + % The actual k-means algorithm. + while any(L ~= L1) + L1 = L; + + % Update clusters + for i = 1:k, l = L==i; + scf = l.*Ix; % signed cluster's frames (1 = pos; 2 = neg) + C(:,i) = (sum(X(:,scf==1),2) - sum(X(:,scf==2),2))/sum(l); + + end + + [M1,Lpos] = min(scoremethod(X,C),[],1); + [M2,Lneg] = min(scoremethod(X,-C),[],1); + [~,Ix] = min([M1;M2]); + L = ones(1,length(Lpos)); L(Ix==1) = Lpos(Ix==1); L(Ix==2) = Lneg(Ix==2); + end + + M = ones(1,length(M1)); M(Ix==1) = M1(Ix==1); M(Ix==2) = M2(Ix==2); + + % Calculate sumD + optVar = sum(M); + + % Is this the best solution so far? If so, remember this one. + if optVar < bestSolution.optVar + bestSolution.optVar = optVar; + bestSolution.L = L; + bestSolution.C = C; + bestSolution.Ix = Ix; + end + end + end + % Return best solution + L = bestSolution.L; + C = bestSolution.C; + Ix = bestSolution.Ix; + + end + + % calculate actual distances + dstpos = distancemethod(X,C); + dstneg = distancemethod(X,-C); + dst = ones(size(dstpos)); dst(:,Ix==1) = dstpos(:,Ix==1); dst(:,Ix==2) = dstneg(:,Ix==2); + + dis = nan(1,k); + sD = nan(1,k); + for i = 1:k + dis(i) = sum(dst(i,L==i)); + sD(i) = std(dst(i,L==i)); + end +end +end + + +%% ------------------------------------------------------------------------ +% DISTANCE FUNCTIONS +% ------------------------------------------------------------------------ +% These functions have been optimised to identify the nearest cluster in a +% vectorised way. +% ------------------------------------------------------------------------ +function dist = sqEuclideanDistance(X,C) + +dist = bsxfun(@plus,dot(X,X,1), sqEuclideanScore(X,C)); + +end + +function dist = CosineDistance(X,C) + +dist = 1 - bsxfun(@rdivide,(X'*C),(sqrt(dot(X,X,1))'*sqrt(dot(C,C,1))))'; + +end + +%% ------------------------------------------------------------------------ +% DISTANCE SCORE FUNCTIONS +% ------------------------------------------------------------------------ +% These functions have been optimised to identify the nearest cluster in a +% vectorised way using the minimum computation possible. Their output +% should be thought of as a score rather than a distance metric. +% ------------------------------------------------------------------------ +function score = sqEuclideanScore(X,C) + +score = bsxfun(@minus, dot(C,C,1), 2*X'*C)'; + +end + +function score = CosineScore(X,C) +% multiply the score by -1 so we optimise for the minimum +% value as opposed to maximum (i.e.: distance, not similarity) + +score = bsxfun(@rdivide, (C'*X)', sqrt(dot(C,C,1)))'*-1; + +end diff --git a/ppicaps/main.m b/ppicaps/main.m new file mode 100644 index 0000000..b3fdd99 --- /dev/null +++ b/ppicaps/main.m @@ -0,0 +1,155 @@ +% ------------------------------------------------------- +% +% main - function to execute PPI-CAPs workflow +% +% Created by: Lorena Freitas +% Last checked: 28.09.2019 +% +% +% ------------------------------------------------------ + +function main + +% ------------------------------------------------------ +% Set up +% ------------------------------------------------------ + +% Distance used for kmeans +dist = 'cosine'; + +% Subjects in analysis +sujets = {'Sujet01', 'Sujet10', 'Sujet16', 'Sujet17', 'Sujet19' ... + 'Sujet21', 'Sujet26', 'Sujet27', 'Sujet31', 'Sujet32', ... + 'Sujet09', 'Sujet14', 'Sujet20', 'Sujet23', 'Sujet25' }; + +ppicapsDir = ''; + +% ------------------------------------------------------ +% Load data for all subjects +% ------------------------------------------------------ + +% Initialise variables +supraFrames=[];subjInfo=[];timeInfo=[];taskInfo=[]; seedSignLabels=[]; + +% Loop through subjects' data to make super subject. The following +% information must have been saved for each subject: the suprathreshold +% frames (selected based on seed's most (de)active moments after whole brain +% deconvolution), the indeces of each suprathreshold frame within the +% original time course (timeInfo), and task labels for each suprathreshold +% frame indicating which task was being performed at that point in time +% (taskInfo); +for i=1:length(sujets) + + % Load subject's variables, paths, etc + thisSubject = initialize_vars(sujets{i}, 'control'); + + % Select suprathresholdFrames + [PPIframes, suprathresholdFramesIx, correspondingTask4Frames] = thresholdFrames(thisSubject, 'PCC', 60); + load([thisSubject.dataDir thisSubject.curSubj '.mat']); % loads PPIframes, suprathresholdFramesIx and correspondingTask4Frames for each subject + + supraFrames = [supraFrames, PPIframes']; % suprathreshold frames from each subject in the shape #voxels x #frames + subjInfo = [subjInfo, ones(1,size( PPIframes',2))*i]; % subject indeces corresponding to each supraframe + timeInfo = [timeInfo, suprathresholdFramesIx]; % each frame's corresponding index within a subject's timecourse + taskInfo = [taskInfo, correspondingTask4Frames]; % task each frame belongs to + + seedSignLabels = [seedSignLabels suprathresholdFramesSeedSign]; + +end + +% Find voxels in and out of Grey Matter +[voxels, nonVoxels] = findGMvoxels; + +% Find matrix rotations for saving nifti files +[voxel_size, voxel_shift] = prepareToSaveNifti; + + +% ------------------------------------------------------ +% Cluster timepoints, using frames as feature vectors +% ------------------------------------------------------ + +% Choose number of clusters (PPI-CAPs) +k = 4; + +% Run kmeans++ on +[PPI_CAPs, PPI_CAPs_Ix] = kmeansCustom(supraFrames(voxels,:), k, 100); + +for i = 1:k + CAPix{i} = find(PPI_CAPs == i); % #frames x 1; index of frames for this CAP + CAPl{i} = PPI_CAPs == i; + scf{i} = CAPl{i}.* PPI_CAPs_Ix; %signed cluster's frames + cluster{i} = supraFrames(:,CAPix{i}); % #voxels x #frames + CAPmean = (sum(supraFrames(:,scf{i}==1),2) - sum(supraFrames(:,scf{i}==2),2))/sum(CAPl{i} ); + CAPmean(nonVoxels,:) = 0; + CAPmean{i} = CAPmean; + + polarityLabels= scf{i}; + polarityLabels = polarityLabels(CAPix{i}); + taskLabels = taskInfo(CAPix{i}); + polarityLabels{i} = polarityLabels; + taskLabels{i} = taskLabels; + seedSigns{i} = seedSignLabels(CAPix{i}); + + + % Temporal normalization, as done by Karahanoglu et al., 2015 + if strcmp(dist, 'cosine') + sd1 = std(cluster{i}, 0, 2); + if (~ismember(0, sd1(voxels)) && ~ismember(Inf, sd1(voxels))) + CAPmean(voxels) = bsxfun(@rdivide, CAPmean(voxels), sd1(voxels)); + CAPmean{i} = CAPmean; + else + warning(['Voxels were not normalised by variance for cap=' num2str(i) ', K=' num2str(k)]); + end + end + + % Spatial nomalization (voxelwise, for visualisation purposes only) + sd2 = std(CAPmean(voxels)); + if (~ismember(0, sd2) && ~ismember(Inf, sd2)) + CAPmean(voxels) = (CAPmean(voxels) - mean(CAPmean(voxels)))/sd2; + CAPmean{i} = CAPmean; + end + + CAP_vol{i} = reshape(CAPmean{i} , 53, 63, 46); % save it back in the original 3D shape + + + % Creates a new NIFTI for the CAP + tmp_cap = make_nii(CAP_vol{i},voxel_size,round(-voxel_shift./voxel_size)); + tmp_cap.hdr.dime.datatype=64; + tmp_cap.hdr.dime.bitpix=64; + + % Saves the cap nifti + capNiftiFile = fullfile(ppicapsDir,['ppicaps_2Level_cap' num2str(i) '_' date '.nii']); + save_nii(tmp_cap,capNiftiFile); + +end + +end + + +% From here on, statistical analysis is performed using the same code as in +% that used in ppicaps_toyExample.m for the statistical analysis part. + + + + +% ------------------------------------------------------ +% HELPER FUNCTIONS +% ------------------------------------------------------ + +% Finds indeces of voxels in or out of Grey Matter (GM) +function [inGM, nonGM] = findGMvoxels + +mask = load([pwd '/dependencies/MNIBrainMask.mat']); % Load MNIBrainMask +nonGM = find(mask.MNIBrainMask(:) < 1); voxels = find(mask.MNIBrainMask(:) == 1); + +end + +function [voxel_size, voxel_shift] = prepareToSaveNifti + +VOI = load([pwd '/dependencies/VOI_PCC_1.mat']); +voxel_size = diag(VOI.xY.spec.mat); +voxel_size = voxel_size(1:end-1)'; + +voxel_shift = VOI.xY.spec.mat(:,4); +voxel_shift = voxel_shift(1:end-1)'; + +end \ No newline at end of file diff --git a/ppicaps/thresholdFrames.m b/ppicaps/thresholdFrames.m new file mode 100644 index 0000000..52fb733 --- /dev/null +++ b/ppicaps/thresholdFrames.m @@ -0,0 +1,138 @@ +% ------------------------------------------------------- +% +% thresholdFrames - function to threshold frames +% +% Created by: Lorena Freitas +% Last checked: 28.09.2019 +% +% +% ------------------------------------------------------ + +function [PPIframes, thresholdedIx, task4frames] = thresholdFrames(thisSubject, seed, percent) + +% ---------------------- +% Variables' set up +% ---------------------- + +% SPM's seed name +VOI_Name = ['VOI_' seed '_1']; + +% Folder where analysis files were stored +b = initialize_vars(thisSubject); +analyze_dir = [b.dataDir 'analysis/']; + +% Total number of frames +nFrames = 179; + +% Load subject's data (deconvData) +load([b.dataDir b.curSubj '_deconvolvedData']); + +% Reshape matrix into 2D to facilitate computations +subjectData2D = reshape(deconvData, [], 340)'; + + +% ---------------------- +% Create seed timecourse +% ---------------------- + +% Obtain seed +VOI = spm_select('FPListRec',analyze_dir,[VOI_Name '.mat$']); VOI = load(VOI); +VOIxyzMNI = VOI.xY.XYZmm; % Dimensions: [3 x #voxels double] +VOIxyzMat = round(VOI.xY.spec.mat \ [VOIxyzMNI; ones(1, length(VOI.xY.XYZmm))]); +VOIxyzMat = VOIxyzMat(1:3,:); % XYZ indeces of seed's voxels in "MATLAB matrix" space +clear VOI_Name VOI VOIxyzMNI + +% Average the seed signal from all voxels +meanSeedSignal = averageSignal(deconvData, VOIxyzMat); + +% ---------------------- +% Create task timecourse +% ---------------------- + +% Repetition time +TR_MW = 2; + +% Onset times for task A (fun scences), in seconds +onset_Fun = [5 22 89 182 288]; + +% Duration of blockes. First row: dur_fun; Second row: dur_sci +dur_MW = [13 8 45 52 23; 4 59 48 54 42]; + +% Get task labels for all frames +[~, taskLabels] = createTaskRegressor(onset_Fun, dur_MW, TR_MW, nScans); + +% Find index of frames to threshold +[thresholdedIx, ~] = thresholdSeed(meanSeedSignal, nFrames, percent); + +% Find suprathreshold frames to keep +PPIframes = subjectData2D(thresholdedIx,:); + +% Find task labels for suprathreshold frames +task4frames = taskLabels(thresholdedIx); + +end + + + + + +% ---------------- +% HELPER FUNCTIONS +% ---------------- + + +function meanSignal = averageSignal(X, VOIxyzMat) + +% Initialise variables +nVoxels = size(VOIxyzMat,2); +signal = zeros(nVoxels, 340); +X = zscore(X); + +% Get signal from voxels +for i=1:nVoxels + x = VOIxyzMat(1,i); y = VOIxyzMat(2,i); z = VOIxyzMat(3,i); + signal(i,:) = squeeze(X(x,y,z,:))'; +end + +% Average signal from voxels +meanSignal = mean(signal); % meanSeedSignal = [1 x # nscan] + +end + + +function [taskRegressor, taskLabels] = createTaskRegressor(onset_Fun, dur_MW, TR_MW, nScans) + +scanTimeFunOnsets = ceil((onset_Fun/TR_MW)); % round up (?) +taskRegressor1 = zeros(1,nscans); + +thisOnset = scanTimeFunOnsets(i); +thisDuration = round(dur_MW(1,i)/TR_MW); +taskRegressor1(thisOnset:(thisOnset+thisDuration)) = 1; + +% Flip zeros and ones +taskRegressor2 = 1 - taskRegressor1; + +% After the 178th scan it's just resting state +taskRegressor2([1,2,175:end]) = 0; +taskRegressor = taskRegressor1; + +taskLabels = taskRegressor; +taskLabels(taskRegressor2==1) = 2; + % create task regressor with 1 / -1 values +taskRegressor(taskRegressor2==1) = -1; + +end + + + +function [thresholdedIx, threshold] = thresholdSeed(meanSeedSignal, nFrames, percent) + +nFrames2keep = floor(nFrames*percent / 100); +[sortedSeedS,sortingIxS] = sort(meanSeedSignal,'descend'); +thresholdedIx = sortingIxS(1:ceil(nFrames2keep/2)); +threshold = min(meanSeedSignal(thresholdedIx)); + +% Now take also moments where the seed is deactive +thresholdedIx = sort([thresholdedIx sortingIxS(end-ceil(nFrames2keep/2):end)]) + +end \ No newline at end of file diff --git a/ppicaps_toyExample.m b/ppicaps_toyExample.m new file mode 100644 index 0000000..76d5c1b --- /dev/null +++ b/ppicaps_toyExample.m @@ -0,0 +1,281 @@ +function ppicaps_toyExample +% +% PPICAPS_TOYEXAMPLE - This function is meant to illustrate how the +% PPI-CAPs time-resolved analysis works. For more information on PPI-CAPs +% please check (Freitas et al., 2019). + +addpath(genpath([pwd '/dependencies'])); + + +% ------------------------------------------------------------------------- +% SET UP SIMULATED EXPERIMENTAL DATA +% ------------------------------------------------------------------------- + +% Let us suppose that we have some experimental task-based data. We would +% like to find out how the relationship between a particular seed activity +% and the organisation of the rest of the brain evolve over time. Are there +% patterns whose correlation with the seed is higher during one specific +% task? And are there others whose relationship with the seed activity is +% constant independently of the task being performed? At what specific +% times do these brain patterns actually happen, throughout the experiment? + +% To answer those questions as per the PPI-CAPs approach, we start by +% selecting time points where the seed is most highly active or deactive. +% This selection can be done based on the percentage of frames you would +% like to keep in the analysis, or based on a threshold of your choice ? as +% we show on the paper, the method is pretty robust to a wide range of +% threshold choises. Since the seed activity must be orthogonal to the task +% timecourse, this will mean that we end up selecting a similar number of +% frames where the seed is active or deactive for each of the tasks. + +% Below, are our true labels for the sign of the seed activity, the task +% timecourse and the ppi (as an element-wise multiplication of the seed and +% task signs), on the suprathreshold frames. That is, we have 24 frames that +% were selected at moments where the seed was highly active or deactive. +seedlabels = [1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1]; +tasklabels = [1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1]; +ppilabels = [1 -1 1 -1 1 -1 -1 1 -1 1 -1 1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 1]; + + +% ------------------------------------------------------------------------- +% SET UP SIMULATED GROUND TRUTH +% ------------------------------------------------------------------------- + +% Let us assume that there are three brain patterns, composed of 5 voxels +% each, as such: +ppicap1 = [-0.7 -0.75 -0.85 -0.9 -0.95]; +ppicap2 = [0.8 0.85 0.9 -0.8 -0.85]; +ppicap3 = [-0.7 0.9 -0.75 0.95 -0.7]; + +% Let's plot them and check them out. +figure; +colors = cbrewer('div', 'Spectral', 50, 'pchip'); caxis([-1.5 1.5]); +imagesc([ppicap1(:),ppicap2(:), ppicap3(:)]);colormap(colors);caxis([-2 2]); +set(gca, 'xtick', 1:3, 'xticklabels',{'Seed', 'Task', 'PPI'},'XTickLabelRotation',90,'FontSize',10); +set(gcf, 'color', 'w', 'pos', [-1562 292 190 473]);colorbar; + + +% At different time points, each of these patterns may be expressed with +% one of two polarities: positive or negative. That is, when ppicap1 is +% expressed with a negative polarity, every voxel has the sign of its +% intensity flipped as compared to its positive polarity. In practical +% terms, this would mean that regions that are active in one polarity are +% de-active in the other polarity, and vice-versa. + +% Moreover, the polarity of each pattern may be random, or it may correlate +% with the activity of a seed region, the timecourse of the task, or an +% interaction of the two, that is: it correlates more with the seed during +% one of the tasks (or it correlates more with the task when the seed is +% either active or deactive). Let's call this interaction a PPI +% (psychophysiological interaction). And we will represent this as shown +% below: + +sp = ppicap1(:); %seedpattern (we will force this to correlate with the seed) +tp = ppicap2(:); %taskpattern (we will force this to correlate with the task) +pp = ppicap3(:); %ppipattern (we will force this to correlate with the ppi) +patternsTimeCourse = [tp, tp, sp, -sp, pp, -pp,... + -tp, -tp, sp, -sp, -pp, pp,... + tp, tp, sp, -sp, pp, -pp,... + -tp, -tp, sp, -sp, -pp, pp]; + +% So let us visualise these patterns. Each one has been labelled with the +% effect (seed / task / ppi, in positive or negative polarity) we have +% forced them to correlate with, for reference purposes. +figure; imagesc(patternsTimeCourse);colormap(colors); +patternsTimeCourseLabels = {'task', 'task', 'seed', '-seed', 'ppi', '-ppi',... + '-task', '-task', 'seed', '-seed', '-ppi', 'ppi',... + 'task', 'task', 'seed', '-seed', 'ppi', '-ppi',... + '-task', '-task', 'seed', '-seed', '-ppi', 'ppi'}; +set(gca, 'xtick', 1:24, 'xticklabels',patternsTimeCourseLabels,'XTickLabelRotation',90,'FontSize',10);set(gcf, 'color', 'w');caxis([-2 2]); + + + +% ------------------------------------------------------------------------- +% STATIONARY ANALYSIS (SiMAP) +% ------------------------------------------------------------------------- + +% This is what a stationary interaction map would look like, when we +% multiply frames by the sign of the corresponding PPI values, and then add +% those together: +figure; imagesc(sum(patternsTimeCourse.*repmat(ppilabels, 5, 1),2));colormap(colors);set(gcf, 'color', 'w');caxis([-16 16]) + +% Note that, in this case, since we only had one pattern that had a PPI +% effect, the stationary analysis was able to recover it perfectly. +% However, if we had two or more different patterns that correlated with +% the PPI time course, the stationary map would return an average of those, +% therefore depicting a pattern that never actually happens in the brain, +% at any point in time!!! This is why we need to proceed with a dynamic +% analysis. + +% ------------------------------------------------------------------------- +% DYNAMIC ANALYSIS ATTEMPTS +% ------------------------------------------------------------------------- + + +% Now let us try to recover the underlying patterns by clustering the +% suprathreshold frames using regular k-means, where k = 3 (since we know +% that this is the exact number of underlying patterns): +[IDX, C] = kmeans(patternsTimeCourse',3,'dist', 'cosine', 'replicates', 100 ); +figure;imagesc(C');colormap(colors); title('K-means (K = 3)');set(gcf, 'color', 'w');caxis([-1 1]); + +% Note that k-means with k = 3 returns patterns that don't even originally +% exist in our dataset, so this would not be a viable option. + + +% Now, let us try again to recover the underlying patterns by clustering +% the frames using k-means with k = 6: +[IDX, C] = kmeans(patternsTimeCourse',6,'dist', 'cosine', 'replicates', 100 ); +figure;imagesc(C');colormap(colors);title('K-means (K = 6)');set(gcf, 'color', 'w'); + +% Note that now we have indeed recover all the different occurences of our +% three patterns. However, when a pattern reappears with a different polarity, +% the two occurences are assigned to different centroids! This would not +% make much sense as the pattern is the same up to a sign flip. More +% importantly, we would need to create and additional metric to +% automatically match centroids that are opposite polarities of the same +% pattern. + + +% ------------------------------------------------------------------------- +% DYNAMIC ANALYSIS WITH PPI-CAPS +% ------------------------------------------------------------------------- + +% So let us instead use PPI-CAPs and cluster using the k-means based on +% modulo-pi cosine distance, to retrieve the flipping of each frame at +% clustering time. +k = 3; +[PPI_CAPs, idx] = kmeanspp(patternsTimeCourse, k, 'cosine', 100); + +for i = 1:k + ppi_cap_vols{i} = patternsTimeCourse(:,find(PPI_CAPs == i)); %instead of frames use test = patternsTimeCourse.*repmat(sign(ppi(idxSupraThresh)), 9,1); + ppi_cap_pols{i} = idx(PPI_CAPs == i); ppi_cap_pols{i}(ppi_cap_pols{i}==2) =-1; + ppi_cap(:,i) = sum(ppi_cap_vols{i}.*repmat(ppi_cap_pols{i}, 5,1),2); +end + +% Now note how we manage to recover the exact patterns we started out with, +% up to a sign flip. Inportantly, the sign itself does not imply that the +% pattern was negative, but that at each given timepoint, each voxel must +% be multiplied by the corresponding sign of the PPI-CAP so we can trully +% know if it is active or deactive at that point in time. +figure;imagesc(ppi_cap); colormap(colors); +title('K-means with modulo-pi dist.');caxis([-16 16]); + + + +% ------------------------------------------------------------------------- +% STATISTICAL ANALYSIS +% ------------------------------------------------------------------------- + +% So now let us suppose we did not know what the group truth was, and let +% us use confusion matrices and permutation testing to check whether each +% PPI-CAP's frames correlate with any of the effects we are trying to +% analyse. + + +seedSigns = seedlabels; +taskSigns = tasklabels; + +maxPerm = 1000; % Maximum number of random permutations +figure; + +% Calculate confusion matrices for each PPI-CAP. +for i = 1:3 + thisCapSeed = seedSigns(PPI_CAPs==i); + thisCapTask = taskSigns(PPI_CAPs==i); + thisCapPPI = thisCapSeed .* thisCapTask; + thisCapClusteringFlip = idx(PPI_CAPs==i); + + % Seed + subplot(3,3,3*i-2); + [seedConfusionMatrix{i}, seed_fkix{i}] = calc_confusionMatrix(thisCapSeed, thisCapClusteringFlip); + imagesc(seedConfusionMatrix{i}); colormap(colors); + + % Task + subplot(3,3,3*i-1); + [taskConfusionMatrix{i}, task_fkix{i}] = calc_confusionMatrix(thisCapTask, thisCapClusteringFlip); + imagesc(taskConfusionMatrix{i});colormap(colors); + + % PPI + subplot(3,3,3*i); + [ppiConfusionMatrix{i}, ppi_fkix{i}] = calc_confusionMatrix(thisCapPPI, thisCapClusteringFlip); + imagesc(ppiConfusionMatrix{i});colormap(colors); + + + % Random Permutations + for nperm = 1: maxPerm + % --- seed + randThisCapSeed = thisCapSeed(randperm(length(thisCapSeed))); + [~, seed_fkix_rnd(i,nperm)]= calc_confusionMatrix(randThisCapSeed, thisCapClusteringFlip); + + % --- task + randThisCapTask = thisCapTask(randperm(length(thisCapTask))); + [~, task_fkix_rnd(i,nperm)]= calc_confusionMatrix(randThisCapTask, thisCapClusteringFlip); + + % --- ppi + randThisCapPPI = thisCapPPI(randperm(length(thisCapPPI))); + [~, ppi_fkix_rnd(i,nperm)]= calc_confusionMatrix(randThisCapPPI, thisCapClusteringFlip); + end + +end + + +% Now let us plot the distributions of fk-indeces obtained through random +% permutations of the (seed; task; and ppi) labels from the frames that +% make up each PPI-CAP. Then, we plot a line for the real data-driven +% fk-index and see where it stands! +figure; +for i = 1:3 + + % Seed + ax2 = subplot(3,3,3*i-2); + h2 = histogram(ax2,seed_fkix_rnd(i,:));title(ax2, 'Null Seed Distribution'); + h2.BinMethod = 'Sturges'; h2.Normalization = 'countdensity'; h2.FaceColor = colors(i,:); h2.NumBins = 5; + line(ax2, [seed_fkix{i}, seed_fkix{i}], get(ax2, 'ylim'), 'LineWidth', 2, 'Color', 'r'); + pSeed = min(((length(find(seed_fkix_rnd(i,:)seed_fkix{i}))+1)/(maxPerm+1))); + text(ax2, 0.7,0.9,['p = ' num2str(pSeed)],'Units','normalized'); + + % Task + ax4 = subplot(3,3,3*i-1); + h2 = histogram(ax4,task_fkix_rnd(i,:));title(ax4, 'Null Task Distribution'); + h2.BinMethod = 'Sturges'; h2.Normalization = 'countdensity'; h2.FaceColor = colors(i,:);h2.NumBins = 5; + line(ax4, [task_fkix{i}, task_fkix{i}], get(ax4, 'ylim'), 'LineWidth', 2, 'Color', 'r'); + pTask = min(((length(find(task_fkix_rnd(i,:)< task_fkix{i}))+1)/(maxPerm+1)),((length(find(task_fkix_rnd(i,:)>task_fkix{i}))+1)/(maxPerm+1))); + text(ax4, 0.7,0.9,['p = ' num2str(pTask)],'Units','normalized'); + + + % PPI + ax3 = subplot(3,3,3*i); + h2 = histogram(ax3, ppi_fkix_rnd(i,:));title(ax3, 'Null PPI Distribution'); + h2.BinMethod = 'Sturges'; h2.Normalization = 'countdensity'; h2.FaceColor = colors(i,:);h2.NumBins = 5; + line(ax3, [ppi_fkix{i}, ppi_fkix{i}], get(ax3, 'ylim'), 'LineWidth', 2, 'Color', 'r'); + pPPI = min(((length(find(ppi_fkix_rnd(i,:)< ppi_fkix{i}))+1)/(maxPerm+1)),((length(find(ppi_fkix_rnd(i,:)>ppi_fkix{i}))+1)/(maxPerm+1))); + text(ax3, 0.7,0.9,['p = ' num2str(pPPI)],'Units','normalized'); + +end + + + +end + + +% This function calculates the confusion matrix for each PPI-CAP. The idea +% is that, if a coactivation pattern has a specific effect (seed / task / +% ppi), the polarity of its composing frames as returned by the clustering +% method should correlate with the sign of that effect. So, for each effect, +% we count the number of frames within that PPI-CAP that are positive when +% the timecourse of that effect is positive, and vice-versa, obtaining a +% confusion matrix. If there is a high correlation between the two, the +% determinant of the matrix will have a high value. The fk-index is thus +% the determinant of the matrix normalised by the total number of frames +% within the PPI-CAPs. +function [confusionMatrix, fk_index] = calc_confusionMatrix(thisCapsEffectLabels, thisCapClusteringFlip) + +confusionMatrix = [sum((thisCapsEffectLabels==1).*(thisCapClusteringFlip ==1)), ... + sum((thisCapsEffectLabels==1).*(thisCapClusteringFlip ==2)); ... + sum((thisCapsEffectLabels==-1).*(thisCapClusteringFlip ==1)), ... + sum((thisCapsEffectLabels==-1).*(thisCapClusteringFlip ==2))]; + +fk_index = det(confusionMatrix)/((numel(thisCapsEffectLabels)/2)^2); + +end +