diff --git a/Utilities/ISFC_thresholding/JOVE_FindExcursions.m b/Utilities/ISFC_thresholding/JOVE_FindExcursions.m new file mode 100644 index 0000000..85dbb54 --- /dev/null +++ b/Utilities/ISFC_thresholding/JOVE_FindExcursions.m @@ -0,0 +1,59 @@ +%% This function computes a distribution of values for the considered +% null data (for each connection pair), and determines (in +% two-tailed fashion) the thresholds past which an excursion of +% connectivity in the data is significant. Then, those moments are returned +% +% Inputs: +% +% - Data and Null1 are the actual and null data (cell arrays n_subject +% with size n_pairs x n_sliding_window_timepoints) acquired during the +% stimulus sessions +% - Null2 contains the purely resting-state recordings, as opposed to the +% resting-state recordings that were acquired within a stimulus-including +% session +% - percentile (in percent) is the level of significance +% +% Outputs: +% +% - Out is a 3D matrix of size n_conn x n_slidingwindow_timepoint x +% subjects, filled with -1 (significant negative ISFC change), 0 (no +% significant change), or +1 (significant positive ISFC change) +function [Out,PCT] = JOVE_FindExcursions(Data,Null,percentile) + + if(size(Data,1) ~= size(Null,1)) + errordlg('The null data and the actual data do not have the same dimensions !!'); + end + + % Number of region pairs + n_pairs = size(Data{1},1); + n_subjects1 = length(Data); + n_subjects2 = length(Null); + + % Aggregates the null (i.e., resting-state) data + Null_marginalized = []; + for s = 1:n_subjects2 + Null_marginalized = [Null_marginalized,Null{s}]; + end + + + % Does the percentile computations + for s = 1:n_subjects1 + PCT{s} = prctile(Null_marginalized,[100-percentile, percentile],2); + end + + % Creates and fills the output, thresholded matrix Out with +-1 at the + % moments when there was a significant ISFC change, for a given subject + % and connection + Out = zeros(n_pairs,size(Data{s},2),n_subjects1); + + for s = 1:n_subjects1 + for i = 1:n_pairs + if sum(Out(i,Data{s}(i,:)>PCT{s}(i,1),s)) > 0 + disp('!!!!!!'); + end + + Out(i,Data{s}(i,:)>PCT{s}(i,1),s) = 1; + Out(i,Data{s}(i,:)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/Utilities/Plotting/cbrewer/cbrewer_preview.jpg b/Utilities/Plotting/cbrewer/cbrewer_preview.jpg new file mode 100644 index 0000000..bd2830a Binary files /dev/null and b/Utilities/Plotting/cbrewer/cbrewer_preview.jpg differ diff --git a/Utilities/Plotting/cbrewer/change_jet.m b/Utilities/Plotting/cbrewer/change_jet.m new file mode 100644 index 0000000..b8d4ecb --- /dev/null +++ b/Utilities/Plotting/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/Utilities/Plotting/cbrewer/colorbrewer.mat b/Utilities/Plotting/cbrewer/colorbrewer.mat new file mode 100644 index 0000000..ec59ef4 Binary files /dev/null and b/Utilities/Plotting/cbrewer/colorbrewer.mat differ diff --git a/Utilities/Plotting/cbrewer/interpolate_cbrewer.m b/Utilities/Plotting/cbrewer/interpolate_cbrewer.m new file mode 100644 index 0000000..e8b5e21 --- /dev/null +++ b/Utilities/Plotting/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/Utilities/Plotting/cbrewer/plot_brewer_cmap.m b/Utilities/Plotting/cbrewer/plot_brewer_cmap.m new file mode 100644 index 0000000..a5cab9e --- /dev/null +++ b/Utilities/Plotting/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/Utilities/Plotting/show_FSmesh.m b/Utilities/Plotting/show_FSmesh.m index da69609..74ce7f4 100755 --- a/Utilities/Plotting/show_FSmesh.m +++ b/Utilities/Plotting/show_FSmesh.m @@ -1,143 +1,143 @@ function myH=show_FSmesh(hemi,meshType,varargin) % SHOW_FSMESH show a freesurfer 3-dimensional hemisphere mesh % % IN: % hemi - a char array containing 'l','r', or 'lr', indicating which hemi- % sphere to display % meshtype - a char array containing mesh type ('white') % 'faceColor': ColorSpec - colour for the faces of the brain mesh % triangles % 'faceAlpha': scalar between 0 and 1 - alpha transparency value to % use for brain mesh triangles (0: fully transparent, 1: fully % opaque) % 'edgeColor': ColorSpec - colour for the edges of the brain mesh % triangles % 'edgeAlpha': scalar between 0 and 1 - alpha transparency value to % use for brain mesh triangle edges (0: fully transparent, 1: fully % opaque) % 'reducePatchVal': scalar, fraction of patches to keep - increasing % from the default value may crash Matlab due to a bug in Matlab. % OUT: % myH: a vector of handles to patch objects % % v1.0 Nov 2009 Jonas Richiardi % - inital release based on read_surf/trisurf snippet by Markus Gschwind % v1.1 2009 Dec 10 Jonas Richiardi % - code clean up % - fix for alpha values crashing in 2-hemispheres display (known bug % 484445): use subdivision algorithm to "downsample" Delaunay triangulation % This is a better solution in this case than those suggested in % http://www.mathworks.com/support/bugreports/484445 % v1.2 2010 Jul Jonas Richiardi % - fix a bug due to reducepatch does not setting CData and FaceVertexCData % after reducing the vertex and face count (despite the doc saying it does % not set the Faces and Vertices properties). Now CData and FaceVertexCData % are in sync. [createFigure,FaceColor,EdgeColor,FaceAlpha,EdgeAlpha,useConvexHull,... reducePatchVal]=process_options(varargin,... 'createFigure',true,'FaceColor',[0.8 0.8 0.8],'EdgeColor',[0.5 0.5 0.5],... 'FaceAlpha',0.2,'EdgeAlpha',0.05,'useConvexHull',false,'reducePatchVal',... 0.015); %% sanity check hemiSupported={'l','r','lr'}; meshTypeSupported={'white'}; if all(cellfun('isempty',strfind(hemiSupported,hemi))) error(['Unsupported hemi type: ' hemi]); end if all(cellfun('isempty',strfind(meshTypeSupported,meshType))) error(['Unsupported mesh type: ' meshType]); end %% setup host-specific options and paths location=whereAmIRunning; if strcmp(location,'jmac') % set root path to FreeSurfer application FSROOTPATH='/Applications/freesurfer/'; else error('Please modify whereAmIRunning.m for your own settings'); end % generate path where to find the surface mesh files FSSURFPATH=fullfile(FSROOTPATH,'subjects','fsaverage','surf'); if exist('read_surf')~=2 addpath(fullfile(FSROOTPATH,'matlab')); end if (createFigure==true) figure; end hold on; %% load surfaces % surfaces are encoded by Delaunay triangulation (region-adjacency graph % of Voronoi tesselation - Delaunay graph is the dual of the Voronoi % diagram (Okabe et al 200)) % faces is the m x 3 triangulation matrix - the set of simplices (triangle) % that make the triangulation. Each row (in simplex-vertex format) encodes % the three vertices a triangular face via indexing into a matrix containing % the 3D space location of vertices % e.g. 1 4 2 says that this triangles is made up of vertices 1, 4, and 2, % whose 3D spatial locations are found in the vertices matrix. myH=zeros(numel(hemi),1); for h=1:numel(hemi) if strcmp(hemi(h),'r') f_surf=['rh.' meshType]; elseif strcmp(hemi(h),'l') f_surf=['lh.' meshType]; else error('Wrong hemi specification. How on earth did you trigger this code path?'); end % read surface using freesurfer function % note we could access matlab functionality for the TriRep class via % trirep(faces+,vertices(:,1)... etc) [vertices,faces]=read_surf(fullfile(FSSURFPATH,f_surf)); vertices=single(vertices); % save memory nVertices = size(vertices,1); if (useConvexHull==true) disp('Computing convex hull...'); myConvHull = convhulln(double(vertices)); myH(h)=trisurf(myConvHull+1,vertices(:,1),vertices(:,2),vertices(:,3),... 'FaceColor',FaceColor,'EdgeColor',EdgeColor); - disp('Doing mesh subdivision...'); + %disp('Doing mesh subdivision...'); %reducepatch(myH(h),0.05); else % plot Delaunay triangulated mesh myH(h)=trisurf(faces+1,vertices(:,1),vertices(:,2),vertices(:,3),... ones(nVertices,1),'FaceColor',FaceColor,'EdgeColor',EdgeColor); % do subdivision to reduce number of faces % another solution would be to set(gca,'DrawMode','Fast'); as per % 484445 bug report - disp('Doing mesh subdivision...'); + %disp('Doing mesh subdivision...'); reducepatch(myH(h),reducePatchVal); % higher may cause crash, still too many faces % CData and FaceVertexCData are now desynchronised, so fix them % manually nf=get(myH(h),'Faces'); nv=get(myH(h),'Vertices'); % 1) fix CData to be 3xsize(Faces,1), all ones or actually FaceColor set(myH(h),'CData',repmat(FaceColor,3,size(nf,1))); % 2) fix FaceVertexCData to be nVertices x 1 set(myH(h),'FaceVertexCData',repmat(FaceColor,size(nv,1),1)); end xlabel('left-right'); ylabel('posterior-anterior'); zlabel('ventral-dorsal'); clear vertices faces; if (h==1) axis equal; view(0,90); end end %set(gca,'DrawMode','Fast'); % set transparency now, with reduced number of faces (does not crash) for h=1:numel(hemi) set(myH(h),'EdgeAlpha',EdgeAlpha,'FaceAlpha',FaceAlpha,... 'DiffuseStrength',1,'SpecularStrength',0); end diff --git a/Utilities/Plotting/show_cm.m b/Utilities/Plotting/show_cm.m index cf92760..fb84532 100755 --- a/Utilities/Plotting/show_cm.m +++ b/Utilities/Plotting/show_cm.m @@ -1,929 +1,929 @@ function [myH,varargout]=show_cm(CM,cod,mode,varargin) % SHOW_CM show a brain connectivity graph with many customisable options % % EXAMPLE: % show_cm(myCM,myCB,3); % % INPUT: % CM: nRegions x nRegions adjacency matrix. Entries indicate graph edge % weights. % cod: a struct for a codebook with fields % .center: 1 x nRegions cell array of 3x1 vectors - MNI coordinates of each region % .id: nRegions x 1 array of integers - code values corresponding to % atlas values (identify each region) % .name: nRegions x 1 cell array of strings - region names % .sname: nRegions x 1 cell array of strings - SHORT region names % .num: a scalar containing the number of regions defined in the % codebook % mode: scalar - % 1 displays adjacency matrix % 2 display 2D projections of graph % 3 shows graph in 3D space % % ** NAME/VALUE PAIRS (incomplete list) % % * OPTIONS TO SELECT WHAT TO PLOT % % 'doPlotConnections': boolean - whether to plot graph edges (true) or % not (default: true) % 'doPlotDegreeSpheres': boolean - whether to plot spheres on vertex % centers representing the degree of the node or a given function (default false) % 'plotDegreeSpheresOnlyAboveThreshold': factor by which degree sphere % must be larger than uniform weight (1/nVertices) in order to be % plotted. 0 - plot all degree spheres (default: 0) / if used together % with FctForColSph, sphere plotted for abs(fct) > threshold % 'nodesXonly': scalar or row vector - show only connections of the specified node(s) % (default: []) % 'pruneTinyConnections': boolean - don't plot connections with very small weight % 'pruningWeight': scalar - weight below which connections should not be % plotted % % % * OPTIONS REGARDING TEXT LABELS % % 'showLabels': logical - show region labels or not (default: true) % 'showWeights': logical - show edge weights or not (defalt: false) % 'myFontSize': scalar between 6 and 50 - size of font for region labels % (default: 12) % 'showTopNlabels': scalar or empty - if scalar, choice of how many % region labels to display - these will be for top-weighted edges. If % empty, label display is either all (showLabels==true) or none. % (default: []) % 'showTopNlabelsBy': whether to take into account edge weights 'edgeW' or vertex % weights 'vertexW'. Default: 'edgeW'; % 'useShortNames': logical - whether to use short names for labels or not % (default: false) % 'txtPosOffset': 1x3 vector - additional label display offset in x,y,z % with respect to the value in center{r}+posOffset (used to make sure % labels can be placed above nodes and edges in certain projections) % (default: 0 0 0) % 'flipLabelsRL': replace 'L' with 'R' at end of region name and % vice-versa (default: false) % % % * OPTIONS TO CONTROL HOW OBJECTS ARE SCALED % % 'linearWeight': scalar - linear scaling for edge weights display % (default: 2) % 'normMax': scalar or empty - if scalar, maximum edge weights with % respect to which all others are scaled. If empty, maximum is % computed from data in CM (default: []) % 'gamma': scalar between 0 and 1 - nonlinear scaling for display of edge % weights (towards 0: small edge weights are inflated, towards 1: small % edge weights are left small) (default: 0.5) % 'f_matD_to_sphS': conversion factor between node degrees and % sphere size for degree plot (default: 1/15) % 'e_matD_to_sphS': exponent for conversion between node degrees and % sphere sizes (default: 2); % 'FctForSizSph': scalar or nVerticesx1, uses same value for all nodes if scalar, % (default: empty) % 'computeDegreesOnFullMat': boolean - compute the weights and degrees % on full (not sparsified through e.g. nodesXonly) adjacency matrix. % Default true % % % * OPTIONS TO CONTROL COLOURS OF VARIOUS ELEMENTS % % 'doAlpha': logical - if true, edges with smaller weights are displayed % more transparent (default: true) % 'doColourCode': logical - if true, edges weight is represented with % colour coding in addition to line thickness (default: true) % 'cmapchoiceE': string - colormap to use for mapping edge weights to % colours. Use standard colormap weights (default: 'autumn') % 'cmapchoiceV': string - colormap to use for mapping Vertex weights to % colours. Use standard colormap weights (default: 'autumn') % 'FctForColSph': nVerticesx1 double, colors spheres on vertices based on % input instead of degree of sphere (default: []) % 'CdataVal': scalar above 0 - value to put in faces of graph edges to % map into a colour map (used if doColourCode==false, useful not to clash % with brain mesh values or other objects) (default: 2) % if = 0, connections colored in gray shade (default: black) % (works only with useSplitcMap = false) % 'GrayShade': scalar [0 1] : colors connections in gray shade (0: black, % 1: white) if 'CdataVal'=0 % 'SphBlack': boolean, if true plots spheres below threshold as gray shade % (c.f. 'GrayShade') instead of omitting them (default: false) % 'BGColour': ColorSpec - background colour for figure (default: 1 1 1) % 'ColAxis': [clim cmax], manually set limits of colorbar (default: []) % 'useSplitCMap': logical - rescale values with respect to splitID so % that different colourmaps can be applied to different subgraphs % ! Experimental code - use freezecolors from Matlab FEX instead % (default: false) % 'splitID': scalar in {1,2} - subgraph identifier, used with % useSplitCMap (default: []) % ! Experimental code % 'alphagamma': scalar between 0 and 1 - further non-linear scaling % applied to edge transparency (allows finer control) (default: .1) % % % * OPTIONS REGARDING 3D CORTICAL MESH % % 'bmesh_show': logical - show brain mesh or not (default: false) % 'bmesh_hemi': string in 'l','r','lr' - select which hemisphere to % display (default: 'lr') % 'bmesh_faceColour': ColorSpec - colour for the faces of the brain mesh % triangles (default: .8 .8 .8) % 'bmesh_faceAlpha': scalar between 0 and 1 - alpha transparency value to % use for brain mesh triangles (0: fully transparent, 1: fully % opaque) (default: .4) % 'bmesh_edgeColour': ColorSpec - colour for the edges of the brain mesh % triangles (default: .5 .5 .5) % 'bmesh_edgeAlpha': scalar between 0 and 1 - alpha transparency value to % use for brain mesh triangle edges (0: fully transparent, 1: fully % opaque) (default: .05) % 'bmesh_useConvexHull': boolean - whether to compute and display the % convex hull rather than the full mesh (default: false) % 'bmesh_convexHullOutlinePlanes': string ('', 'a','s','c'): if not empty, % draw an outline of the convex hull on the corresponding plane (axial, % sagittal, or coronal). Needs the bmesh_useConvexHull option to work. % (default: '') % 'bmesh_convexHullOutlinePlanesLineWidth': width of convex hull outline % plot (default: 1) % 'bmesh_reducePatchVal': scalar - fraction of patches to retain for % drawing brain meshes - increasing this from default of 0.015 may crash % Matlab. % % % * GENERAL AND MISC OPTIONS % % 'createFigure: logical - whether a new figure should be created or not % (default: true) % 'posOffset': 3x1 vector - display offset in x,y,z with respect to the % value in center{r} (default: 0 0 0) % 'Title': add a title to the plot (default: []) % % % OUTPUT: % myH: a cell array of graphics handles. % myH{1}: scalar, contains handle to the figure if createFigure is % true, empty otherwise % myH{2}: vector, contains 1 or 2 handles for the patch objects % representing the hemispheric meshes or convex hull thereof % if bmesh_hemi is set, otherwise is empty % myH{3}: vector, contains 1 to 3 handles for the plot3 line objects % representing the a,s,c projections of the convex hull of the % hemispheric meshes, if bmesh_convexHullOutlinePlanes is set, % otherwise is empty % myH{4}: vector, contains handles to all surface objects (degree % spheres and connections) % myH{5}: vector, contains handles to all text objects % varargout(1): % wCM: a reweighted and thresholded CM corresponding to % what show_CM displays. the diagonal corresponds to region % sphere sizes. % % VERSION HISTORY: % v1.0 (date unknown) Dimitri Van de Ville % - initial release % v2.0 oct-nov 2009 Jonas Richiardi % - allow for adjacency matrix plotting (zero values = no edge) % - nicer color mapping % - 3D view with patches and transparency control % - all weights rescaled to max 1 % - optional node label and weights display % - gamma exponent for non-linearity selectable % - alpha transparency on/off % - single-colour or colour mapped to link strength % - returns figure handles % - many more... % v2.0.1 2009 Nov 10 Jonas Richiardi % - code cleanup for combined mesh + connection view % v2.1 2010 April 21 Jonas Richiardi % - convex hull computation and plotting, and outline of convex hull. % v3.0 jul 2010 Nora Leonardi % - print only connections of 1 node % v3.1 jul 2010 Jonas Richiardi % - cleaner handle return structure. Separate handles for figure, hemispheric % meshes, hull outline projectsion, surface-plotted objects, and text, to % allow the use of separate renderers (painters / opengl) when exporting % to eps. This allows the generate of high-quality text and surface % objects for publications, and keeps editorial offices and compositors % happy (see jHighQRenderingToEPS.m). % - other minor cleanups and improvements % v3.1.1 jul 2010 Jonas Richiardi % - bugfix on 3.1.0: handles of only half coronal and axial outlines were % available % v3.1.2 Jul 2010 Jonas Richiardi % - don't floor and display very small weight connections, instead don't % display them at all % v 3.2 Mar 2011 JR % - minor bugfixes % v3.2.1 Nov 2011 JR % - partial merge with Nora's branch - FctForSizSph % v3.2.2 Nov 2011 NL % - FctForSizSph can now be used to control sphere size in parallel % with FctForColSph controlling the colour % v3.2.3 Mar 2013 JR % - fix logical condition code for nodesXonly % v3.3 April 2013 JR % - reorganised options, now sorted alphabetically % - harmonised capitalisation between incoming options and internal names % - improved documentation, but still work in progress % v3.4 Aug 2013 NL % - degrees calculated using abs(CM) % - removed Cdataval and Grayshade warning as called with default inputs % - changed pp.FctForSizSph: input is scalar or nVerticesx1 % showWeights off by default, text() is a slow function [pp.alphagamma,pp.BGColour,pp.bmesh_convexHullOutlinePlanes,... pp.bmesh_convexHullOutlinePlanesLineWidth,pp.bmesh_edgeAlpha,... pp.bmesh_edgeColour,pp.bmesh_faceAlpha,pp.bmesh_faceColour,... pp.bmesh_hemi,pp.bmesh_reducePatchVal, pp.bmesh_show,... pp.bmesh_useConvexHull,pp.CdataVal,pp.cmapchoiceE,pp.cmapchoiceV,... pp.ColAxis,pp.computeDegreesOnFullMat,pp.createFigure,pp.doAlpha, ... pp.doColourCode,pp.doPlotConnections,pp.doPlotDegreeSpheres,... pp.e_matD_to_sphS, pp.f_matD_to_sphS,pp.FctForColSph,... pp.FctForSizSph,pp.flipLabelsRL,pp.gamma,pp.GrayShade,... pp.linearWeight,pp.myFontSize,pp.nodesXonly,pp.normMax,... pp.plotDegreeSpheresOnlyAboveThreshold,pp.posOffset,... pp.pruneTinyConnections,pp.pruningWeight,pp.showLabels,... pp.showTopNlabels,pp.showTopNlabelsBy,pp.showWeights,pp.SphBlack,... pp.splitID,pp.Title,pp.txtPosOffset,pp.useShortNames,... pp.useSplitCMap]=... process_options(varargin,'alphagamma',0.1,'BGColour',[1 1 1],... 'bmesh_convexHullOutlinePlanes','',... 'bmesh_convexHullOutlinePlanesLineWidth',1,'bmesh_edgeAlpha',0.05,... 'bmesh_edgeColour',[0.5 0.5 0.5],'bmesh_faceAlpha',0.05,... 'bmesh_faceColour',[0.8 0.8 0.8],'bmesh_hemi','lr',... 'bmesh_reducePatchVal',0.015,'bmesh_show',false,... 'bmesh_useConvexHull',false,'CdataVal',2,'cmapchoiceE','autumn',... 'cmapchoiceV','autumn','ColAxis',[],'computeDegreesOnFullMat',true,... 'createFigure',true,'doAlpha',true,'doColourCode',true,... 'doPlotConnections',true,'doPlotDegreeSpheres',false,... 'e_matD_to_sphS',2,'f_matD_to_sphS', 1/15,'FctForColSph',[],... 'FctForSizSph',[],'flipLabelsRL',false,'gamma',0.5,... 'GrayShade',0,'linearWeight',2,'myFontSize',12,'nodesXonly',[],... 'normMax',[],'plotDegreeSpheresOnlyAboveThreshold',0,... 'posOffset',[0 0 0]','pruneTinyConnections',false,... 'pruningWeight',0.01,'showLabels',true,... 'showTopNlabels',[],'showTopNlabelsBy','edgeW','showWeights',false,... 'SphBlack',false,'splitID',[],'Title',[],'txtPosOffset',[0 0 0]',... 'useShortNames',false,'useSplitCMap',false); %% CONSTANTS % handle cell array identification constants in myH H_F=1; % handle to the figure H_H=2; % handle to hemispheric meshes H_O=3; % handle to brain outlines H_S=4; % handle to surface objects H_T=5; % handle to text objects %% sanity check % check codebook allCBfields=fieldnames(cod); needCBfields={'id','num','name','sname','center'}; for f=1:numel(needCBfields) if all(cellfun('isempty',strfind(allCBfields,needCBfields{f}))) error(['Codebook is lacking field: ' needCBfields{f}]); end end % check params are within permissible ranges if (pp.gamma <0) || (pp.gamma > 1) error('gamma must be between 0 and 1'); end if (pp.linearWeight<=0) error('linearWeight must be above 0'); end if (pp.bmesh_faceAlpha<0) || (pp.bmesh_faceAlpha>1) error('bmesh_faceAlpha must be between 0 and 1'); end if (pp.myFontSize<7) || (pp.myFontSize>50) error('bmesh_faceAlpha must be between 7 and 50'); end if ~isempty(pp.showTopNlabels) && pp.showLabels && ... ((pp.showTopNlabels <1) || (pp.showTopNlabels > cod.num)) error('showTopNlabels must be between 1 and cod.num'); end if pp.CdataVal<0 error('CdataVal must be greater than 0'); end if ~isempty(pp.bmesh_convexHullOutlinePlanes) if pp.bmesh_useConvexHull==false || pp.bmesh_show==false error('bmesh_useConvexHull and bmesh_show must be true to use bmesh_convexHullOutlinePlanes'); end for p=1:numel(pp.bmesh_convexHullOutlinePlanes) if all(cellfun('isempty',strfind({'a','s','c'},pp.bmesh_convexHullOutlinePlanes(p)))) error('bmesh_convexHullOutlinePlanes must be empty, a, s, or c'); end end end % if pp.FctForSizSph && isempty(pp.FctForColSph) % error('FctForColSph cannot be empty if FctForSizSph is true'); % end % check nodes to plot exist if ~isempty(pp.nodesXonly) if any(pp.nodesXonly > cod.num | pp.nodesXonly==0) error('elements of nodesXonly must be valid indices of CM (1 to %d)', cod.num); end end % check splitmap if (pp.useSplitCMap==true) if isempty(pp.splitID) error('Must specify splitID when using useSplitCMap'); elseif ~ismember(pp.splitID,[1 2]); error('splitID must be 1 or 2'); end end %if (pp.alphagamma <0) || (pp.alphagamma > 1) % error('alphagamma must be between 0 and 1'); %end % check offsets if ~all(size(pp.posOffset)==[3 1]) if ~all(size(pp.posOffset)==[1 3]) error('Please supply column vectors for posOffset'); else pp.posOffset=pp.posOffset'; end end if ~all(size(pp.txtPosOffset)==[3 1]) if ~all(size(pp.txtPosOffset)==[1 3]) error('Please supply column vectors for txtPosOffset'); else pp.txtPosOffset=pp.txtPosOffset'; end end % check CDataVal=0 if GrayShade not empty % if pp.CdataVal~=0 && ~isempty(pp.GrayShade) % warning('GrayShade not used since CdataVal > 0'); % end if pp.SphBlack==true && isempty(pp.plotDegreeSpheresOnlyAboveThreshold) error('Please specifiy a threshold below which spheres are plotted in gray shade: plotDegreeSpheresOnlyAboveThreshold'); end % compute complementary colour for text pp.compColour=abs([1 1 1]-pp.BGColour); % replace long names by short ones if needed if (pp.useShortNames==true) cod.name=cod.sname; end if ~isempty(pp.FctForSizSph) && length(pp.FctForSizSph)==1 pp.FctForSizSph=pp.FctForSizSph*ones(size(CM,1),1); end %% cut adjacency matrix if only connetions of 1 node to be shown if ~isempty(pp.nodesXonly) A=zeros(cod.num); Nnodes=length(pp.nodesXonly); for k=1:Nnodes colx=CM(:,pp.nodesXonly(k)); rowx=CM(pp.nodesXonly(k),:); % could also reuse colx as CM sym A(:,pp.nodesXonly(k))=colx; A(pp.nodesXonly(k),:)=rowx; end CMold = CM; % keep old matrix to calc node degrees CM = A; else CMold = CM; end %% main selector % display adjacency matrix if mode==1, if (pp.createFigure==true) figure; caxis([-1 1]); end imagesc(CM); h=get(gcf,'CurrentAxes'); if (pp.showLabels==true) set(h,'Ytick',[1:length(CM)]); ml=0; for iter=1:cod.num; ml=max(ml,length(cod.name{iter})); end; dispOffset=5; ml=ml+dispOffset; % save space for staggering % pre-fill str with blanks str=repmat([' '],[cod.num ml]); for iter=1:cod.num; if mod(iter,2)==0 offset=dispOffset; else offset=1; end str(iter,offset:offset+length(cod.name{iter})-1)=cod.name{iter}; end; set(h,'YTickLabel',str,'FontSize',8); end caxis([-1 1]); axis square; axis tight; % plot connections in 2D or 3D else % normalise by max if isempty(pp.normMax) maxVal=max(max(CM)); maxValOld=max(max(CMold)); else maxVal=pp.normMax; maxValOld=pp.normMax; end CM=CM./maxVal; CMold=CMold./maxValOld; % get column vectors for r=1:cod.num cod.center{r}=cod.center{r}(:); end % 2D brain view if mode==2, warning('returned handle vector is not compatible with show_cm v3.1'); myH=zeros(3,1); % display myH(1)=brain_plot(CM,cod,[1 3],pp); title('coronal'); xlabel('left-right'); ylabel('ventral-dorsal'); myH(2)=brain_plot(CM,cod,[1 2],pp); title('axial'); xlabel('left-right'); ylabel('posterior-anterior'); myH(3)=brain_plot(CM,cod,[2 3],pp); title('sagittal'); xlabel('posterior-anterior'); ylabel('ventral-dorsal'); % 3D brain view elseif mode==3 if (pp.bmesh_show==true) nHemis=numel(pp.bmesh_hemi); else nHemis=0; end myH=cell(5,1); if nHemis==0 myH{H_H}=[]; else myH{H_H}=zeros(nHemis,1); % hemispheric meshes or hulls end if (pp.createFigure==true) myH{H_F}=figure; else myH{H_F}=[]; end hold on; caxis([0 1]); % show brain mesh if (pp.bmesh_show==true) myH{H_H}=show_FSmesh(pp.bmesh_hemi,... 'white','createFigure',false,'FaceColor',pp.bmesh_faceColour,... 'FaceAlpha',pp.bmesh_faceAlpha,'EdgeColor',pp.bmesh_edgeColour,... 'EdgeAlpha',pp.bmesh_edgeAlpha,'useConvexHull',pp.bmesh_useConvexHull, ... 'reducePatchVal',pp.bmesh_reducePatchVal); % setup mesh colours for hi=1:nHemis set(myH{H_H}(hi),'FaceColor',pp.bmesh_faceColour,... 'FaceAlpha',pp.bmesh_faceAlpha,'EdgeColor',pp.bmesh_edgeColour,... 'EdgeAlpha',pp.bmesh_edgeAlpha); end end % plot graph [myH{H_S},myH{H_T},wCM]=brain_plot3d(CM,cod,pp,CMold); varargout(1)={wCM}; myH{H_T}(end+1)=xlabel('left-right','Color',pp.compColour); myH{H_T}(end+1)=ylabel('posterior-anterior','Color',pp.compColour); myH{H_T}(end+1)=zlabel('ventral-dorsal','Color',pp.compColour); if ~isempty(pp.Title), myH{H_T}(end+1)=title(pp.Title); end % plot hull outline planes if needed if (pp.bmesh_useConvexHull==true) && ~isempty(pp.bmesh_convexHullOutlinePlanes) myH{H_O}=[]; % zeros(numel(pp.bmesh_convexHullOutlinePlanes)*nHemis,1); Xs=cell(nHemis,1); Ys=cell(nHemis,1); Zs=cell(nHemis,1); kas=cell(nHemis,1);kss=cell(nHemis,1);kcs=cell(nHemis,1); for hi=1:nHemis % get all 3D vertices locations in convex hull xd = get(myH{H_H}(hi),'XData'); yd = get(myH{H_H}(hi),'YData'); zd = get(myH{H_H}(hi),'ZData'); Xs{hi}=[xd(1,:) xd(2,:) xd(3,:)]; Ys{hi}=[yd(1,:) yd(2,:) yd(3,:)]; Zs{hi}=[zd(1,:) zd(2,:) zd(3,:)]; % now we can compute the convex hull of the projections kas{hi}=convhull(Xs{hi},Ys{hi}); % convex hull of the axial projection kss{hi}=convhull(Ys{hi},Zs{hi}); % sagittal kcs{hi}=convhull(Xs{hi},Zs{hi}); % coronal % we could also create slice plane with values - see "exploring volumes with slice %planes" end %figure; scatter3(Xs{1},Ys{1},Zs{1},'k'); hold on; %scatter3(Xs{2},Ys{2},Zs{2},'k'); axis equal; sagHasBeenPlotted=false; for hi=1:nHemis % delete old 3D convex hull graphics object and set its % handle to nana delete(myH{H_H}(hi)); myH{H_H}(hi)=nan; hold on; % and now plot this stuff to check for p=1:numel(pp.bmesh_convexHullOutlinePlanes) switch pp.bmesh_convexHullOutlinePlanes(p) case 'a' myH{H_O}(end+1)=plot3(Xs{hi}(kas{hi}),Ys{hi}(kas{hi}),zeros(numel(kas{hi}),1),... 'k--','LineWidth',pp.bmesh_convexHullOutlinePlanesLineWidth); % plot axial case 's' if ~sagHasBeenPlotted % avoid plotting twice - hemispheres one behind the other myH{H_O}(end+1)=plot3(zeros(numel(kss{hi}),1),Ys{hi}(kss{hi}),Zs{hi}(kss{hi}),'k--',... 'LineWidth',pp.bmesh_convexHullOutlinePlanesLineWidth); % plot sagittal sagHasBeenPlotted=true; end case 'c' myH{H_O}(end+1)=plot3(Xs{hi}(kcs{hi}),zeros(numel(kcs{hi}),1),Zs{hi}(kcs{hi}),'k--',... 'LineWidth',pp.bmesh_convexHullOutlinePlanesLineWidth); % plot coronal otherwise error('bmesh_convexHullOutlinePlanes must be a,s,c. There is no way you should be here!'); end end % all convex hull outline planes end % all hemispheres end % plot convex hull outline plane else error('Unknown visualisation mode'); end end function myH=brain_plot(CM,cod,coor,pp) % 2D warning('show_cm:ancientCode','Ancient code. Check myH return format and plotting details.'); myH=figure;hold on; for iter=1:cod.num, scatter(cod.center{iter}(coor(1)),cod.center{iter}(coor(2)),'o','MarkerEdgeColor',pp.compColour); if (pp.showLabels==true) text(cod.center{iter}(coor(1)),cod.center{iter}(coor(2)),... cod.name{iter},'Interpreter','none','HorizontalAlignment',... 'center','Color',pp.compColour); end end; for iter=1:cod.num, for iter2=iter+1:cod.num, if ~isnan(CM(iter,iter2)) && CM(iter,iter2)~=0 % compute display params myW=abs(CM(iter,iter2)); edgeWeight=pp.linearWeight*myW^pp.gamma; %colVal=exp(2*(1-myEdges(e).Weight)); colVal=1-myW; if colVal < 0.2, colVal=0.2; end if colVal > 0.9, colVal=0.9; end h=plot([cod.center{iter}(coor(1)) cod.center{iter2}(coor(1))],[cod.center{iter}(coor(2)) cod.center{iter2}(coor(2))],'b-'); set(h,'LineWidth',edgeWeight); set(h,'Color',[colVal colVal colVal]); if (pp.showWeights==true) midPt=mean([cod.center{iter}([coor(1) coor(2)]) ... cod.center{iter2}([coor(1) coor(2)])],2); text(midPt(1),midPt(2),num2str(edgeWeight,'%1.2f')); end end; end; end; axis equal; hold off; function [myHs,myHt,wCM]=brain_plot3d(CM,cod,pp,CMold) % IN: % CM: square matrix, a (potentially sparsified) adjacency matrix for a graph % cod: struct, a codebook annotating regions % pp: struct, plotting parameters % CMold: square matrix, original full adjacency matrix % % OUT: % myHs: vector, handles to all surface objects % myHt: vector, handles to all text objects myHs=[]; % surface objects myHt=[]; % text objects nVertices=size(CM,1); wCM=zeros(nVertices)*nan; %posOffset=0; %100; nCylPts=20; % if (pp.createFigure==true) % mH=figure; % hold on; % caxis([0 1]); % else % mH=[]; % hold on; % caxis([0 1]); % end set(gca,'Color',pp.BGColour,'XColor',pp.compColour,'YColor',pp.compColour,... 'ZColor',pp.compColour); set(gcf,'Color',pp.BGColour); %% plot region centroids and vertex degrees if pp.computeDegreesOnFullMat==true % check symmetry if (all(all(CMold==CMold'))) CMsym=CMold; else CMsym=(CMold+CMold')/2; end else if (all(all(CM==CM'))) CMsym=CM; else CMsym=(CM+CM')/2; end end % compute soft degrees %myDegrees=sum(triu(CMsym,1),2); myDegrees=sum(abs(CMsym)); % use abs weight [sXbase,sYbase,sZbase] = sphere(30); % coord for sphere -> draw sphere with surf (filled) or mesh (grid) for iter=1:cod.num, cx=cod.center{iter}(1)+pp.posOffset(1); cy=cod.center{iter}(2)+pp.posOffset(2); cz=cod.center{iter}(3)+pp.posOffset(3); if isempty(pp.FctForSizSph) thisDeg=myDegrees(iter)+1; else % use value i fct to size thisDeg=pp.FctForSizSph(iter)+1; end %plot3(cx,cy,cz,'o','MarkerEdgeColor',pp.compColour./3,... % 'MarkerSize',10); % thisDeg if pp.doPlotDegreeSpheres % XXX the else/if below seem unnecessarily complicated if abs(thisDeg-1)>=pp.plotDegreeSpheresOnlyAboveThreshold*(1/nVertices) foo_plotDS=true; else foo_plotDS=false; end if ~isempty(pp.FctForColSph) && pp.SphBlack==true % OVERRIDE ABOVE, plots all foo_plotDS=true; end % if isempty(pp.FctForColSph) % use degree to decide if to plot % if abs(thisDeg-1)>=pp.plotDegreeSpheresOnlyAboveThreshold*(1/nVertices) % foo_plotDS=true; % else % foo_plotDS=false; % end % elseif ~isempty(pp.FctForColSph) % if pp.SphBlack==true % plots all % foo_plotDS=true; % else % uses degree as above % if ((thisDeg-1)>=pp.plotDegreeSpheresOnlyAboveThreshold*(1/nVertices)) % foo_plotDS=true; % else % foo_plotDS=false; % end % %if (abs(pp.FctForColSph(iter))>=pp.plotDegreeSpheresOnlyAboveThreshold) % % foo_plotDS=true; % %else % % foo_plotDS=false; % %end % end % end % else % foo_plotDS=true; % end if foo_plotDS==true thisDegS=(thisDeg*pp.f_matD_to_sphS)^pp.e_matD_to_sphS;%^4; wCM(iter,iter)=thisDegS; sX=(thisDegS*sXbase)+cx; sY=(thisDegS*sYbase)+cy; sZ=(thisDegS*sZbase)+cz; if ~isempty(pp.FctForColSph) % plot function on vertex instead of degree myCData=pp.FctForColSph(iter); if ~isempty(pp.ColAxis), caxis(pp.ColAxis); end else myCData=((thisDegS/pp.f_matD_to_sphS)/(nVertices+1)); % normalise between 0 and 1 : max degree is nVertices end if pp.SphBlack==true && ~(abs(pp.FctForColSph(iter))>=pp.plotDegreeSpheresOnlyAboveThreshold) myCoData=ones([size(sZ,1) size(sZ,2) 3]); myCoData=myCoData*pp.GrayShade; % same color as connections else myCoData=repmat(myCData,size(sZ)); end myHs(end+1)=surf(sX,sY,sZ,'FaceAlpha',1.0,'CData',myCoData,'edgecolor','none'); end end end; %% colormap % some colormaps need to be inverted so that low edge weights are bright % and high vals are dark invCMs={'autumn','bone','hot','gray','copper','pink'}; if all(cellfun('isempty',strfind(invCMs,pp.cmapchoiceV))) colormap(pp.cmapchoiceV); else colormap(flipud(colormap(pp.cmapchoiceV))); end if strfind(pp.cmapchoiceV,'hsv') colormap([hsv;hsv]); end freezeColors; %% plot connections and weights edgeWeights=zeros(size(CM),'single'); plotThisConnection=false; -disp('Processing nodes...'); +%disp('Processing nodes...'); for iter=1:cod.num, if (mod(iter,10)==0) - fprintf('%s ',num2str(iter)); + %fprintf('%s ',num2str(iter)); %disp(['Processing node ' num2str(iter) '...']); end for iter2=iter+1:cod.num, if ~isnan(CM(iter,iter2)) && CM(iter,iter2)~=0 % compute display params myWraw=abs(CM(iter,iter2)); myW=myWraw^pp.gamma; edgeWeight=pp.linearWeight*myW; %^2; edgeWeights(iter,iter2)=edgeWeight; %colVal=myW; colVal = CM(iter,iter2) ^ pp.gamma; if pp.pruneTinyConnections==true && myW 0.99, colVal=sign(colVal)*0.99; end if pp.doAlpha==true alphaVal=myWraw^pp.alphagamma; % ^3 for more rapid decrease if alphaVal0.99, alphaVal=0.99; end else alphaVal=0.99; end %h=plot3([cod.center{iter}(1) cod.center{iter2}(1)],[cod.center{iter}(2) cod.center{iter2}(2)],[cod.center{iter}(3) cod.center{iter2}(3)],'b-'); %set(h,'LineWidth',edgeWeight*3); %set(h,'Color',[colVal colVal colVal]); [fooX, fooY, fooZ] = cylinder2P(edgeWeight, ... nCylPts,[cod.center{iter}]'+pp.posOffset',... [cod.center{iter2}]'+pp.posOffset'); if (pp.doColourCode==true) if pp.useSplitCMap==true % map split 1 to 0-0.5 and split 2 to 0.5-1 colVal=(colVal/2)+0.5*(pp.splitID-1); % what falls on 0.5 is mapped down for split 1, up for % split 2 if (colVal==0.5) if (pp.splitID==1) colVal=colVal-eps; else colVal=colVal+eps; end end end myCdata=repmat(colVal,2,nCylPts); else % assign a fixed scalar CDATA value to this vertex if pp.useSplitCMap==true % assign different CDATA values to different splits myCdata=repmat(pp.CdataVal+1*(pp.splitID-1),2,nCylPts); else % all vertices use same CDATAvalue (pp.CdataVal) if pp.CdataVal==0 % show connections as black myCdata=ones([2 nCylPts 3]); myCdata=myCdata*pp.GrayShade; else myCdata=repmat(pp.CdataVal,2,nCylPts); end end end if (pp.doPlotConnections==true) % each side is defined by XData YData ZData (all 2xnCylPts), where % 3D points are defined by X(i,j) Y(i,j) Z(i,j) % and lines joining them are X(:,j) Y(:,j) Z(:,j) % sides can be drawn by surface(X(:,1:2),Y(:,1:2),Z(:,1:2)) % force double for cData to avoid complaint - maybe % input data is single... myHs(end+1)=surf(fooX,fooY,fooZ,'FaceAlpha',alphaVal,'EdgeAlpha',alphaVal','CData',double(myCdata)); else %myHs(end+1)=[]; end if (pp.showWeights==true) midPt=mean([cod.center{iter}+pp.posOffset cod.center{iter2}+pp.posOffset],2); myHt(end+1)=text(midPt(1),midPt(2),midPt(3),[num2str(colVal,'%1.2f')]); % ... % '=>' num2str(edgeWeight,'%1.2f')]); end end % plotThisConnection end; end; end; %% plot labels % compute subset to show (long-winded way) if ~isempty(pp.showTopNlabels) && pp.showLabels labelsToShow=zeros(cod.num,1); if strcmp(pp.showTopNlabelsBy,'edgeW') edgeWeights_v=jUpperTriMatToVec(edgeWeights,1); edgeWeights_v_srt=sort(edgeWeights_v,'descend'); cutoffW=edgeWeights_v_srt(pp.showTopNlabels); if cutoffW==0 % there might be fewer non-zero edges than showTopNlabels warning('Not enough non-zero weights for pp.showTopNlabels parameter, selecting fewer'); nonZeroWeights=find(edgeWeights_v_srt); if ~isempty(nonZeroWeights) cutoffW=edgeWeights_v_srt(nonZeroWeights(end)); else cutoffW=Inf; end end % point out edges whose nodes are worth labelling connsToShow=edgeWeights>=cutoffW; % ugly for i=1:cod.num for j=i:cod.num if (connsToShow(i,j)==true) labelsToShow(i)=1; labelsToShow(j)=1; end end end elseif strcmp(pp.showTopNlabelsBy,'vertexW') % XXX hack myDegrees_srt=sort(myDegrees,'descend'); % myDegrees cutoffW=myDegrees_srt(pp.showTopNlabels); labelsToShow=myDegrees >= cutoffW; else error('Unknown option for plotting labels by'); end else labelsToShow=ones(cod.num,1); end % do the display if (pp.showLabels==true) for iter=1:cod.num if (labelsToShow(iter)==1) processedLabel=cod.name{iter}; if pp.flipLabelsRL if strcmp(processedLabel(end),'L') processedLabel(end)='R'; else processedLabel(end)='L'; end end myHt(end+1)=text(double(cod.center{iter}(1)+pp.posOffset(1)+pp.txtPosOffset(1)),... double(cod.center{iter}(2)+pp.posOffset(2)+pp.txtPosOffset(2)),... double(cod.center{iter}(3)+pp.posOffset(3)+pp.txtPosOffset(3)),... processedLabel,... 'FontWeight','normal','FontSize',pp.myFontSize,... 'Interpreter','none','HorizontalAlignment','center', ... 'Color',pp.compColour); end end end -fprintf('%s\n','Done'); +%fprintf('%s\n','Done'); axis equal; hold off; %view(66,10); grid off; lighting flat; shading flat; % colormap edge if all(cellfun('isempty',strfind(invCMs,pp.cmapchoiceE))) colormap(pp.cmapchoiceE); else colormap(flipud(colormap(pp.cmapchoiceE))); end diff --git a/Utilities/Preprocessing/RemoveNaNRegions.m b/Utilities/Preprocessing/RemoveNaNRegions.m index 2dd14f5..469095a 100644 --- a/Utilities/Preprocessing/RemoveNaNRegions.m +++ b/Utilities/Preprocessing/RemoveNaNRegions.m @@ -1,47 +1,43 @@ %% This function verifies whether a given atlas time course is acceptable % or not % % Inputs: % % - TCin is the input cell array containing the data of each subject in a % cell (of size n_ROI x n_timepoints) % % Outputs: % % - TCout is the same array, having removed the bad regional time courses % - idx_toremove summarizes which regional indices have been removed % % Written by Thomas Bolton, last checked on July 24th function [TCout,idx_toremove] = RemoveNaNRegions(TCin) % Indices of the ROIs where in at least one subject, there are NaN % values idx_toremove = []; - n_set = length(TCin); - n_ROI = size(TCin{1},1); + n_ROI = size(TCin,1); % If there is at least a NaN value in a given time course, we include % it in the ones to remove - for i = 1:n_set - if ~isempty(TCin{i}) - for j = 1:n_ROI - if(sum(isnan(TCin{i}(j,:))) > 0) - idx_toremove = [idx_toremove,j]; - end + if ~isempty(TCin) + for j = 1:n_ROI + if(sum(isnan(TCin(j,:))) > 0) + idx_toremove = [idx_toremove,j]; end end end % Makes each ROI appear only once in the vector idx_toremove = unique(idx_toremove); % Only keeps the time courses that were deemed OK by the check - for i = 1:n_set - if ~isempty(TCin{i}) - TCout{i} = TCin{i}(~ismember(1:n_ROI,idx_toremove),:); - else - TCout{i} = NaN; - end + if ~isempty(TCin) + TCout = TCin(~ismember(1:n_ROI,idx_toremove),:); + else + TCout = NaN; end + end \ No newline at end of file