diff --git a/Plotting/Create_CAP_colorbar.m b/Plotting/Create_CAP_colorbar.m new file mode 100755 index 0000000..be85c38 --- /dev/null +++ b/Plotting/Create_CAP_colorbar.m @@ -0,0 +1,56 @@ +%% Creates a nice looking colorbar for the CAP display +function [Ah,h] = Create_CAP_colorbar(absmin,absmax,absstep,colT,lab,Ah,Ori,CB1,CB2,n_steps) + + H_range = [absmin absmax]; % The colormap is symmetric around zero + + % Set the Min/Max T-values for alpha coding + A_range = [0 1]; + + % Set the labels for the colorbar + hue_label = lab; + + colrange = linspace(absmin,absmax,256); + + switch Ori + case 'Horizontal' + + y = linspace(A_range(1), A_range(2), 256); + % x represents the range in alpha (abs(t-stats)) + x = linspace(H_range(1), H_range(2), 256); + % y represents the range in hue (beta weight difference) + [X,Y] = meshgrid(x,y); % Transform into a 2D matrix + + h=imagesc(x,y,X,'Parent',Ah); + axis(Ah,'xy'); % Plot the colorbar + set(Ah, 'Xcolor', 'k', 'Ycolor', 'k','YTickLabel','','YTick',[],'XTick',absmin:absstep:absmax,'FontSize',8); + set(Ah, 'XAxisLocation', 'bottom'); + xlabel(Ah,hue_label,'FontSize',8); + + A = ones(size(X)); + A(abs(X) < colT) = 0; + A = reshape(A,256,256); + + case 'Vertical' + + x = linspace(A_range(1), A_range(2), 256); + % x represents the range in alpha (abs(t-stats)) + y = linspace(H_range(1), H_range(2), 256); + % y represents the range in hue (beta weight difference) + [X,Y] = meshgrid(x,y); % Transform into a 2D matrix + + h=imagesc(x,y,Y,'Parent',Ah); + axis(Ah,'xy'); % Plot the colorbar + set(Ah, 'Xcolor', 'k', 'Ycolor', 'k','XTickLabel','','XTick',[],'YTick',absmin:absstep:absmax,'FontSize',8); + set(Ah, 'YAxisLocation', 'right'); + ylabel(Ah,hue_label,'FontSize',8); + + A = ones(size(Y)); + A(abs(Y) < colT) = 0; + A = reshape(A,256,256); + end + + tmp_cmap = cbrewer(CB1,CB2,n_steps); + colormap(Ah,flipud(tmp_cmap)); + + set(h,'AlphaData',A); +end \ No newline at end of file diff --git a/Plotting/MakeViolin.m b/Plotting/MakeViolin.m new file mode 100755 index 0000000..6dc59ee --- /dev/null +++ b/Plotting/MakeViolin.m @@ -0,0 +1,56 @@ +%% Makes a violin plot with overlapped boxplot, tuning colors +% Y and O are vectors containing the data to plot for the two conditions +% Color has light colors for background (first cell) and dark colors for +% average (second cell), for maximum 4 populations +function [box,h_viol,ah] = MakeViolin(Y,ah,Lab,YLabel,Color,n_pop,n_states) + + % Range of the plot + Max = max(max((Y))); + Min = min(min((Y))); + Max = Max + 0.15*max(abs(Min),Max); + Min = Min - 0.15*max(abs(Min),Max); + + % Plots the distribution + ylim(ah,[Min,Max]); + + % Colors for the background + Col_final = num2cell(Color{1},2); + Col_final = Col_final(1:n_pop); + Col_final = repmat(Col_final,n_states,1); + + % Colors for the foreground + Col_final2 = num2cell(Color{2},2); + Col_final2 = Col_final2(1:n_pop); + Col_final2 = repmat(Col_final2,n_states,1); + + Lab_final = {}; + + for i = 1:n_states + Lab_final = [Lab_final, repmat({Lab{i}},1,n_pop)]; + end + + h_viol = distributionPlot(ah,Y','showMM',0,'color',Col_final,'yLabel',YLabel,'xNames',Lab_final); + h_viol=0; + hold(ah,'on'); + box = boxplot(Y','Parent',ah,'Labels',Lab_final); + + h_box = findobj(box,'Tag','Box'); + set(h_box,'color','k','LineWidth',2); + + h_median = findobj(box,'Tag','Median'); + h_outliers = findobj(box,'Tag','Outliers'); + + for i = 1:n_pop + for j = 1:n_states + + idx_oi = i+(j-1)*n_pop; + + set(h_median(idx_oi),'color',Col_final2{idx_oi},'LineWidth',2); + set(h_outliers(idx_oi),'Marker','o','MarkerFaceColor',Col_final2{n_pop*n_states-idx_oi+1},... + 'MarkerEdgeColor',Col_final2{idx_oi},'MarkerSize',3); + + end + end + + set(ah,'Box','off'); +end \ No newline at end of file diff --git a/Plotting/cbrewer/cbrewer.m b/Plotting/cbrewer/cbrewer.m new file mode 100644 index 0000000..26be891 --- /dev/null +++ b/Plotting/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/Plotting/cbrewer/cbrewer_preview.jpg b/Plotting/cbrewer/cbrewer_preview.jpg new file mode 100644 index 0000000..bd2830a Binary files /dev/null and b/Plotting/cbrewer/cbrewer_preview.jpg differ diff --git a/Plotting/cbrewer/change_jet.m b/Plotting/cbrewer/change_jet.m new file mode 100644 index 0000000..b8d4ecb --- /dev/null +++ b/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/Plotting/cbrewer/colorbrewer.mat b/Plotting/cbrewer/colorbrewer.mat new file mode 100644 index 0000000..ec59ef4 Binary files /dev/null and b/Plotting/cbrewer/colorbrewer.mat differ diff --git a/Plotting/cbrewer/interpolate_cbrewer.m b/Plotting/cbrewer/interpolate_cbrewer.m new file mode 100644 index 0000000..e8b5e21 --- /dev/null +++ b/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/Plotting/cbrewer/plot_brewer_cmap.m b/Plotting/cbrewer/plot_brewer_cmap.m new file mode 100644 index 0000000..a5cab9e --- /dev/null +++ b/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/Plotting/distributionPlot.m b/Plotting/distributionPlot.m new file mode 100755 index 0000000..352cc69 --- /dev/null +++ b/Plotting/distributionPlot.m @@ -0,0 +1,902 @@ +function handles = distributionPlot(varargin) +%DISTRIBUTIONPLOT creates violin plots for convenient visualization of multiple distributions +% +% SYNOPSIS: handles = distributionPlot(data,propertyName,propertyValue,...) +% handles = distributionPlot(ah,...) +% +% INPUT data : m-by-nData array of values, or vector of grouped data (use +% the 'groups' property to specify the grouping variable), or +% cell array of length nData. +% The cell array can either contain vectors with values, or +% m-by-2 arrays with [bins,counts] if you want to determine the +% histograms by yourself (m can be different between cell +% elements). Note that arrays inside cells with any +% other shape than m-by-2 are reshaped to vector an a warning is +% thrown (DISTRIBUTIONPLOT:AUTORESHAPE). +% +% DISTRIBUTIONPLOT accepts the following propertyName/propertyValue +% pairs (all are optional): +% +% distWidth : width of distributions; ideally between 0 and 1. +% 1 means that adjacent distributions might touch. Default: 0.9 +% variableWidth : If true, the width of the distribution changes, +% reflecting the shape of the histogram of the data. If false, +% the distribution is only encoded by color levels. Default: true +% color : uniform coloring of histograms. Supply either a color +% string ('r'), or a truecolor vector ([1 0 0]). Use a +% cell array of length nData to specify one color per +% distribution. Default: 'k' +% If variableWidth is set to false, a colormap is generated that +% goes from white to the chose color (or from black, if +% invert==true). +% If both 'color', and 'colormap' are specified, 'colormap' takes +% precedence. +% colormap : colormap used to describe the distribution (first row +% corresponds to bins with least data, last row corresponds to +% bins with most data (invert the grayscale colormap to have +% black indicate the most data). +% Supply a cell array of length nData to color distributions +% individually. Note that using multiple colormaps means that +% the colorbar doesn't contain much useful information. +% Default: [] +% Colormap will index into the figure colormap, which will be +% modified by distributionPlot. This is done to allow editing the +% distributions in e.g. Adobe Illustrator. +% If both 'color', and 'colormap' are specified, 'colormap' takes +% precedence. +% globalNorm : normalization for bin width (x-direction) +% 0 : every histogram is normalized individually so that the +% maximum bin width is equal to distWidth. This is best +% suited to comparing distribution shapes. Default. +% 1 : histograms are normalized such that equal bin width +% reports equal numbers of counts per bin. +% 2 : histograms are normalized so that the relative areas +% covered by the histograms reflect the relative total number +% of data points. +% 3 : histograms areas are normalized so that relative densities +% are the same across histograms. Thus, if +% data = {rand(100,1),rand(500,1)}, +% then +% distributionPlot(data,'globalNorm',2,'histOpt',0,'divFactor',10) +% shows the left histogram 5x as wide as the right, while +% distributionPlot(data,'globalNorm',3,'histOpt',0,'divFactor',10) +% displays both histograms equally wide, since each bin +% contains ~10% of the data. +% Options 1 and 2 produce similar results if the bins are spaced +% equally for the distributions. Options 0 and 3 produce similar +% results if the data are drawn from the same distributions. +% Note that colormaps currently always report the number of data +% points per bin; 'globalNorm' only applies to the distribution +% shape. +% +% groups : grouping variable for grouped data. Grouping will be +% resolved by calling grp2idx, and unless xNames have +% been supplied, group names determine the x-labels. +% If the grouping variable is numeric, group labels also +% determine x-values, unless the parameter xValues has +% been specified. +% histOpt : histogram type to plot +% 0 : use hist command (no smoothing, fixed number of +% bins) +% 1 : smoothened histogram using ksdensity with +% Normal kernel. Default. +% 1.1: smoothened histogram using ksdensity where the +% kernel is robustly estimated via histogram.m. +% Normal kernel. +% 2 : histogram command (no smoothing, automatic +% determination of thickness (y-direction) of bins) +% divFactor : Parameter dependent on histOpt. If... +% histOpt == 0: divFactor = # of bins. Default: 25. +% Alternatively, pass a vector which will be +% interpreted as bin centers. +% histOpt == 1: divFactor decides by how much the default +% kernel-width is multiplied in order to avoid an +% overly smooth histogram. Default: 1/2 +% histOpt == 2: divFactor decides by how much the +% automatic bin width is multiplied in order to have +% more (<1) or less (>1) detail. Default: 1 +% addSpread : if 1, data points are plotted with plotSpread. +% distWidth is ideally set to 0.95 +% This option is not available if the data is supplied as +% histograms. +% Please download plotSpread.m separately from the File +% Exchange using the link in the remarks +% showMM : if 1, mean and median are shown as red crosses and +% green squares, respectively. This is the default +% 2: only mean +% 3: only median +% 4: mean +/- standard error of the mean (no median) +% 5: mean +/- standard deviation (no median) +% 6: draw lines at the 25,50,75 percentiles (no mean) +% 0: plot neither mean nor median +% xValues: x-coordinate where the data should be plotted. +% If xValues are given, "distWidth" is scaled by the median +% difference between adjacent (sorted) x-values. Note that +% this may lead to overlapping distributions. Default: +% 1:nData +% xNames : cell array of length nData containing x-tick names +% (instead of the default '1,2,3') +% xMode : if 'auto', x-ticks are spaced automatically. If 'manual', +% there is a tick for each distribution. If xNames is +% provided as input, xMode is forced to 'manual'. Default: +% 'manual'. +% NOTE: SPECIFYING XNAMES OR XVALUES OR XMODE WILL ERASE PREVIOUS +% LABELS IF PLOTTING INTO EXISTING AXES +% yLabel : string with label for y-axis. Default : '' +% If empty and data is histograms, ylabel is set to 'counts' +% invert : if 1, axes color is changed to black, and colormap is +% inverted. +% histOri: Orientation of histogram. Either 'center', 'left', or +% 'right'. With 'left' or 'right', the left or right half of +% the standard violin plot is shown. Has no effect if +% variableWidth is false. Default: center +% xyOri : orientation of axes. Either 'normal' (=default), or +% 'flipped'. If 'flipped', the x-and y-axes are switched, so +% that violin plots are horizontal. Consequently, +% axes-specific properties, such as 'yLabel' are applied to +% the other axis. +% widthDiv : 1-by-2 array with [numberOfDivisions,currentDivision] +% widthDiv allows cutting the stripe dedicated to a single +% distribution into multible bands, which can be filled with +% sequential calls to distributionPlot. This is one way +% to compare two (or more) sequences of distributions. See +% example below. +% ah : axes handle to plot the distributions. Default: gca +% +% OUTPUT handles : 1-by-4 cell array with patch-handles for the +% distributions, plot handles for mean/median, the +% axes handle, and the plotSpread-points handle +% +% +% EXAMPLES +% %--Distributions contain more information than boxplot can capture +% r = rand(1000,1); +% rn = randn(1000,1)*0.38+0.5; +% rn2 = [randn(500,1)*0.1+0.27;randn(500,1)*0.1+0.73]; +% rn2=min(rn2,1);rn2=max(rn2,0); +% figure +% ah(1)=subplot(3,4,1:2); +% boxplot([r,rn,rn2]) +% ah(2)=subplot(3,4,3:4); +% distributionPlot([r,rn,rn2],'histOpt',2); % histOpt=2 works better for uniform distributions than the default +% set(ah,'ylim',[-1 2]) +% +% %--- additional options +% +% data = [randn(100,1);randn(50,1)+4;randn(25,1)+8]; +% subplot(3,4,5) +% +% %--- defaults +% distributionPlot(data); +% subplot(3,4,6) +% +% %--- show density via custom colormap only, show mean/std, +% distributionPlot(data,'colormap',copper,'showMM',5,'variableWidth',false) +% subplot(3,4,7:8) +% +% %--- auto-binwidth depends on # of datapoints; for small n, plotting the data is useful +% % note that this option requires the additional installation +% % of plotSpread from the File Exchange (link below) +% distributionPlot({data(1:5:end),repmat(data,2,1)},'addSpread',true,'showMM',false,'histOpt',2) +% +% %--- show quantiles +% subplot(3,4,9),distributionPlot(randn(100,1),'showMM',6) +% +% %--- horizontal orientation +% subplot(3,4,10:11), +% distributionPlot({chi2rnd(3,1000,1),chi2rnd(5,1000,1)},'xyOri','flipped','histOri','right','showMM',0), +% xlim([-3 13]) +% +% %--- compare distributions side-by-side (see also example below) +% % plotting into specified axes will throw a warning that you can +% % turn off using " warning off DISTRIBUTIONPLOT:ERASINGLABELS " +% ah = subplot(3,4,12); +% subplot(3,4,12),distributionPlot(chi2rnd(3,1000,1),'histOri','right','color','r','widthDiv',[2 2],'showMM',0) +% subplot(3,4,12),distributionPlot(chi2rnd(5,1000,1),'histOri','left','color','b','widthDiv',[2 1],'showMM',0) +% +% %--Use globalNorm to generate meaningful colorbar +% data = {randn(100,1),randn(500,1)}; +% figure +% distributionPlot(data,'globalNorm',true,'colormap',1-gray(64),'histOpt',0,'divFactor',[-5:0.5:5]) +% colorbar +% +% %--Use widthDiv to compare two series of distributions +% data1 = randn(500,5); +% data2 = bsxfun(@plus,randn(500,5),0:0.1:0.4); +% figure +% distributionPlot(data1,'widthDiv',[2 1],'histOri','left','color','b','showMM',4) +% distributionPlot(gca,data2,'widthDiv',[2 2],'histOri','right','color','k','showMM',4) +% +% %--Christmas trees! +% x=meshgrid(1:10,1:10); +% xx = tril(x); +% xx = xx(xx>0); +% figure +% hh=distributionPlot({xx,xx,xx},'color','g','addSpread',1,'histOpt',2,'showMM',0); +% set(hh{4}{1},'color','r','marker','o') +% END +% +% REMARKS To show distributions as clouds of points (~beeswarm plot), +% and/or to use the option "addSpread", please download the +% additional function plotSpread.m from the File Exchange +% http://www.mathworks.com/matlabcentral/fileexchange/37105-plot-spread-points-beeswarm-plot +% +% I used to run ksdensity with the Epanechnikov kernel. However, +% for integer data, the shape of the kernel can produce peaks +% between the integers, which is not ideal (use histOpt=2 for +% integer valued data). +% +% A previous iteration of distributionPlot used the input +% specifications below. They still work to ensure backward +% compatibility, but are no longer supported or updated. +% handles = distributionPlot(data,distWidth,showMM,xNames,histOpt,divFactor,invert,addSpread,globalNorm) +% where distWidth of 1 means that the maxima +% of two adjacent distributions might touch. Negative numbers +% indicate that the distributions should have constant width, i.e +% the density is only expressed through greylevels. +% Values between 1 and 2 are like values between 0 and 1, except +% that densities are not expressed via graylevels. Default: 1.9 +% +% +% SEE ALSO histogram, ksdensity, plotSpread, boxplot, grp2idx +% + +% created with MATLAB ver.: 7.6.0.324 (R2008a) on Windows_NT +% +% created by: Jonas Dorn; jonas.dorn@gmail.com +% DATE: 08-Jul-2008 +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%==================================== +%% TEST INPUT +%==================================== + +% set defaults +def.xNames = []; +def.showMM = 1; +def.distWidth = 0.9; +def.histOpt = 1; +def.divFactor = [25,2,1]; +def.invert = false; +def.colormap = []; +def.color = 'k'; +def.addSpread = false; +def.globalNorm = false; +def.variableWidth = true; +def.groups = []; +def.yLabel = ''; +def.xValues = ''; +def.xMode = 'manual'; +def.histOri = 'center'; +def.xyOri = 'normal'; +def.widthDiv = [1 1]; +isHistogram = false; %# this parameter is not set by input + + +if nargin == 0 || isempty(varargin{1}) + error('not enough input arguments') +end + +% check for axes handle +if ~iscell(varargin{1}) && isscalar(varargin{1}) == 1 && ... + ishandle(varargin{1}) && strcmp(get(varargin{1},'Type'),'axes') + ah = varargin{1}; + data = varargin{2}; + varargin(1:2) = []; + newAx = false; + + +else + ah = gca; + data = varargin{1}; + varargin(1) = []; + newAx = true; +end + +% check for current axes limits. Set NaN if the axes have no children +% yet - we need that in case we're building a complicated set of +% distributions +if ~isempty(get(ah,'children')) + xAxLim = xlim; + yAxLim = ylim; +else + [xAxLim,yAxLim] = deal([NaN NaN]); +end + +fh = get(ah,'Parent'); + +% check data. If not cell, convert +if ~iscell(data) + [nPoints,nData] = size(data); + data = mat2cell(data,nPoints,ones(nData,1)); +else + % get nData + data = data(:); + nData = length(data); + % make sure all are vectors + badCol = ~cellfun(@isvector,data) & ~cellfun(@isempty,data); + if any(badCol) + nCols = cellfun(@(x)(size(x,2)),data(badCol)); + if all(nCols==2) + % bins,counts + isHistogram = true; + else + warning('DISTRIBUTIONPLOT:AUTORESHAPE',... + 'Elements %s of the cell array are not vectors. They will be reshaped automatically',... + num2str(find(badCol)')); + data(badCol) = cellfun(@(x)(x(:)),data(badCol),'UniformOutput',false); + end + end +end + +parserObj = inputParser; +parserObj.FunctionName = 'distributionPlot'; +stdWidth = 1; % scaling parameter for variableWidth with uneven x-values +% check whether we're dealing with pN/pV or straight arguments +if ~isempty(varargin) && ~ischar(varargin{1}) && ~isstruct(varargin{1}) + % use old format + % distWidth,showMM,xNames,histOpt,divFactor,invert,addSpread,globalNorm + def.distWidth = 1.9; + parserObj.addOptional('distWidth',def.distWidth); + parserObj.addOptional('showMM',def.showMM); + parserObj.addOptional('xNames',def.xNames); + parserObj.addOptional('histOpt',def.histOpt); + parserObj.addOptional('divFactor',def.divFactor); + parserObj.addOptional('invert',def.invert); + parserObj.addOptional('addSpread',def.addSpread); + parserObj.addOptional('globalNorm',def.globalNorm); + parserObj.addOptional('groups',def.groups); + parserObj.addOptional('yLabel',def.yLabel); + parserObj.addOptional('color',def.color); + + + parserObj.parse(varargin{:}); + opt = parserObj.Results; + % fill in defaults that are not supported in the old version of the + % code + opt.colormap = []; + opt.variableWidth = true; + opt.histOri = 'center'; + opt.xValues = []; + opt.xMode = 'auto'; + opt.xyOri = 'normal'; + opt.widthDiv = [1 1]; + + % overwrite empties with defaults - inputParser considers empty to be a + % valid input. + fnList = fieldnames(opt); + for fn = fnList' + if isempty(opt.(fn{1})) + opt.(fn{1}) = def.(fn{1}); + end + end + + + % fix a few parameters + if opt.distWidth > 1 + opt.distWidth = opt.distWidth - 1; + else + opt.colormap = 1-gray(128); + end + if opt.distWidth < 0 + opt.variableWidth = false; + opt.distWidth = abs(opt.distWidth); + end + + if ~isempty(opt.xNames) + opt.xMode = 'manual'; + end + + +else + defNames = fieldnames(def); + for dn = defNames(:)' + parserObj.addParamValue(dn{1},def.(dn{1})); + end + + + parserObj.parse(varargin{:}); + opt = parserObj.Results; + + % if groups: deal with data + if ~isempty(opt.groups) + [idx,labels,vals] = grp2idx(opt.groups); + % convert data to cell array + data = accumarray(idx,data{1},[],@(x){x}); + nData = length(data); + % if not otherwise provided, use group labels for xnames + if isempty(opt.xNames) + opt.xNames = labels; + if ~iscell(opt.xNames) + opt.xNames = num2cell(opt.xNames); + end + end + if isnumeric(vals) && isempty(opt.xValues) + opt.xValues = vals; + end + + end + + if ~ischar(opt.xyOri) || ~any(ismember(opt.xyOri,{'normal','flipped'})) + error('option xyOri must be either ''normal'' or ''flipped'' (is ''%s'')',opt.xyOri); + end + + + + +end +% common checks + +% default x-values: 1:n +if isempty(opt.xValues) + opt.xValues = 1:nData; +elseif length(opt.xValues) ~= nData + error('please supply as many x-data values as there are data entries') +elseif length(opt.xValues) > 1 % only check for scale if more than 1 value + % scale width + stdWidth = median(diff(sort(opt.xValues))); + opt.distWidth = opt.distWidth * stdWidth; +end + +if ~isscalar(opt.divFactor) && length(opt.divFactor) == 3 && all(opt.divFactor==def.divFactor) + opt.divFactor = opt.divFactor(floor(opt.histOpt)+1); +end +if isHistogram + opt.histOpt = 99; + if isempty(opt.yLabel) + opt.yLabel = 'counts'; + end +end + + + +% check colors/colormaps: do we need to expand colormap? +if ~iscell(opt.colormap) + opt.colormap = {opt.colormap}; +end +if ~iscell(opt.color) + opt.color = {opt.color}; +end +for iColor = 1:length(opt.color) + if ischar(opt.color{iColor}) + opt.color{iColor} = colorCode2rgb(opt.color{iColor}); + end +end + +% expand - if only single colormap specified, we expand only once +if ~opt.variableWidth + missingColormaps = find(cellfun(@isempty,opt.colormap)); + for iMissing = missingColormaps(:)' + + endColor = opt.color{max(iMissing,length(opt.color))}; + % normally, we go from white to color + cmap = zeros(128,3); + for rgb = 1:3 + cmap(:,rgb) = linspace(1,endColor(rgb),128); + end + opt.colormap{iMissing} = cmap; + + end +end + +% if we have colormaps, we need to create a master which we add to the +% figure. Invert if necessary, and expand the cell array to nData +colormapLength = cellfun(@(x)size(x,1),opt.colormap); +if any(colormapLength>0) + + colormap = cat(1,opt.colormap{:}); + if opt.invert + colormap = 1-colormap; + end + set(fh,'Colormap',colormap) + if length(opt.colormap) == 1 + opt.colormap = repmat(opt.colormap,nData,1); + colormapLength = repmat(colormapLength,nData,1); + colormapOffset = zeros(nData,1); + singleMap = true; + else + colormapOffset = [0;cumsum(colormapLength(1:end-1))]; + singleMap = false; + end + +else + + colormapLength = zeros(nData,1); + if length(opt.color) == 1 + opt.color = repmat(opt.color,nData,1); + end + if opt.invert + opt.color = cellfun(@(x)1-x,opt.color,'uniformOutput',false); + end +end + + +% set hold on +holdState = get(ah,'NextPlot'); +set(ah,'NextPlot','add'); + +% if new axes: invert +if newAx && opt.invert + set(ah,'Color','k') +end + +%=================================== + + + +%=================================== +%% PLOT DISTRIBUTIONS +%=================================== + +% assign output +hh = NaN(nData,1); +[m,md,sem,sd] = deal(nan(nData,1)); +if opt.showMM == 6 + md = nan(nData,3,3); % md/q1/q3, third dim is y/xmin/xmax +end + +% get base x-array +% widthDiv is a 1-by-2 array with +% #ofDivs, whichDiv +% The full width (distWidth) is split into +% #ofDivs; whichDiv says which "stripe" is active +xWidth = opt.distWidth/opt.widthDiv(1); +xMin = -opt.distWidth/2; +xLow = xMin + xWidth * (opt.widthDiv(2)-1); +xBase = [-xWidth;xWidth;xWidth;-xWidth]/2; +xOffset = xLow + xWidth/2; + +% b/c of global norm: loop twice +plotData = cell(nData,2); + +% loop through data. Prepare patch input, then draw patch into gca +for iData = 1:nData + currentData = data{iData}; + % only plot if there is some finite data + if ~isempty(currentData(:)) && any(isfinite(currentData(:))) + + switch floor(opt.histOpt) + case 0 + % use hist + [xHist,yHist] = hist(currentData,opt.divFactor); + + case 1 + % use ksdensity + + if opt.histOpt == 1.1 + % use histogram to estimate kernel + [dummy,x] = histogram(currentData); %#ok + if length(x) == 1 + % only one value. Make fixed distribution + dx = 0.1; + yHist = x; + xHist = sum(isfinite(currentData)); + else + dx = x(2) - x(1); + + % make sure we sample frequently enough + x = min(x)-dx:dx/3:max(x)+dx; + [xHist,yHist] = ksdensity(currentData,x,'kernel','normal','width',dx/(1.5*opt.divFactor)); + end + else + + % x,y are switched relative to normal histogram + [xHist,yHist,u] = ksdensity(currentData,'kernel','normal'); + % take smaller kernel to avoid over-smoothing + if opt.divFactor ~= 1 + [xHist,yHist] = ksdensity(currentData,'kernel','normal','width',u/opt.divFactor); + end + end + + % modify histogram such that the sum of bins (not the + % integral under the curve!) equals the total number of + % observations, in order to be comparable to hist + xHist = xHist/sum(xHist)*sum(isfinite(currentData)); + + case 2 + % use histogram - bar heights are counts as in hist + [xHist,yHist] = histogram(currentData,opt.divFactor,0); + case 99 + % bins,counts already supplied + xHist = currentData(:,2)'; + yHist = currentData(:,1)'; + end + plotData{iData,1} = xHist; + plotData{iData,2} = yHist; + end +end + +goodData = find(~cellfun(@isempty,plotData(:,1))); +% get norm +switch opt.globalNorm + case 3 + % #3 normalizes relative densities + xNorm(goodData) = cellfun(@(x)min(diff(x)),plotData(goodData,2)); + xNorm(goodData) = xNorm(goodData) .* cellfun(@sum,plotData(goodData,1))'; + maxNorm(goodData) = cellfun(@max,plotData(goodData,1)); + xNorm(goodData) = xNorm(goodData)*max(maxNorm(goodData)./xNorm(goodData)); + + case 2 + % #2 should normalize so that the integral of the + % different histograms (i.e. area covered) scale with the + % respective sum of counts across all bins. Requires evenly spaced + % histograms at the moment + xNorm(goodData) = cellfun(@(x)min(diff(x)),plotData(goodData,2)); + maxNorm(goodData) = cellfun(@max,plotData(goodData,1)); + xNorm(goodData) = xNorm(goodData)*max(maxNorm(goodData)./xNorm(goodData)); + case 1 + xNorm(goodData) = max(cat(2,plotData{:,1})); + case 0 + xNorm(goodData) = cellfun(@max,plotData(goodData,1)); +end + + +for iData = goodData' + + % find current data again + currentData = data{iData}; + + xHist = plotData{iData,1}; + yHist = plotData{iData,2}; + + % find y-step + dy = min(diff(yHist)); + if isempty(dy) + dy = 0.1; + end + + % create x,y arrays + nPoints = length(xHist); + xArray = repmat(xBase,1,nPoints); + yArray = repmat([-0.5;-0.5;0.5;0.5],1,nPoints); + + + % x is iData +/- almost 0.5, multiplied with the height of the + % histogram + if opt.variableWidth + + + tmp = xArray.*repmat(xHist,4,1)./xNorm(iData); + + switch opt.histOri + case 'center' + % we can simply use xArray + xArray = tmp; + case 'right' + % shift everything to the left + delta = tmp(1,:) - xArray(1,:); + xArray = bsxfun(@minus,tmp,delta); + case 'left' + % shift everything to the right + delta = tmp(1,:) - xArray(1,:); + xArray = bsxfun(@plus,tmp,delta); + end + + xArray = xArray + opt.xValues(iData); + + else + xArray = xArray + iData; + end + + % add offset (in case we have multiple widthDiv) + xArray = xArray + xOffset; + + + % yData is simply the bin locations + yArray = repmat(yHist,4,1) + dy*yArray; + + % add patch + vertices = [xArray(:),yArray(:)]; + faces = reshape(1:numel(yArray),4,[])'; + + if colormapLength(iData) == 0 + colorOpt = {'FaceColor',opt.color{iData}}; + else + % calculate index into colormap + if singleMap + % use scaled mapping so that colorbar is meaningful + if opt.globalNorm > 0 + colorOpt = {'FaceVertexCData',xHist','CDataMapping','scaled','FaceColor','flat'}; + else + colorOpt = {'FaceVertexCData',xHist'/xNorm(iData),'CDataMapping','scaled','FaceColor','flat'}; + end + + else + idx = round((xHist/xNorm(iData))*(colormapLength(iData)-1))+1; + colorOpt = {'FaceVertexCData',idx'+colormapOffset(iData),'CDataMapping','direct','FaceColor','flat'}; + end + end + + + switch opt.xyOri + case 'normal' + hh(iData)= patch('Vertices',vertices,'Faces',faces,'Parent',ah,colorOpt{:},'EdgeColor','none'); + case 'flipped' + hh(iData)= patch('Vertices',vertices(:,[2,1]),'Faces',faces,'Parent',ah,colorOpt{:},'EdgeColor','none'); + end + + if opt.showMM > 0 + if isHistogram + [m(iData),sem(iData)] = weightedStats(currentData(:,1),currentData(:,2),'w'); + sd(iData) = sem(iData) * sqrt(sum(currentData(:,2))); + % weighted median: where we're at middle weight + % may need some tweaking + goodCurrentData = sortrows(currentData(all(isfinite(currentData),2),:),1); + weightList = cumsum(goodCurrentData(:,2)); + weightList = weightList / weightList(end); + md(iData) = goodCurrentData(find(weightList>0.5,1,'first'),1); + else + m(iData) = nanmean(currentData); + md(iData) = nanmedian(currentData); + sd(iData) = nanstd(currentData); + sem(iData) = sd(iData)/sqrt(sum(isfinite(currentData))); + end + + if opt.showMM == 6 + % read quantiles - "y"-value, plus x-start-stop + % re-use md array which allows using a loop below instead of + % lots of copy-paste + % md array is md/q1/q3, with third dimension y/xmin/xmax + + md(iData,2,1) = prctile(currentData,25); + md(iData,3,1) = prctile(currentData,75); + + for qq = 1:3 + % find corresponding y-bin + yLoc = repmat(... + any(yArray>md(iData,qq,1),1) & any(yArray<=md(iData,qq,1),1),... + [4 1]); + % look up corresponding x-values. Note that there is a bit + % of a risk that the line will be exactly between two very + % different bins - but if we make the line longer, it will + % be ugly almost all the time + md(iData,qq,2) = min( xArray( yLoc ) ); + md(iData,qq,3) = max( xArray( yLoc ) ); + end + + end + end +end % loop + +sh = []; +if opt.addSpread + if isHistogram + disp('Option addSpread is unavailable if data is supplied as histograms. Call plotSpread separately') + else + % add spread + try + sh = plotSpread(ah,data,'xValues',opt.xValues,'xyOri',opt.xyOri); + set(sh{1},'color',[0,128,255]/255); + catch me + if strcmp(me.identifier,'MATLAB:UndefinedFunction') + error('plotSpread not found. Please download it from the Matlab File Exchange') + else + rethrow(me) + end + end + end +end + +mh = [];mdh=[]; +if opt.showMM + % plot mean, median. Mean is filled red circle, median is green square + % I don't know of a very clever way to flip xy and keep everything + % readable, thus it'll be copy-paste + switch opt.xyOri + case 'normal' + if any(opt.showMM==[1,2]) + mh = plot(ah,opt.xValues+xOffset,m,'+r','Color','r','MarkerSize',12); + end + if any(opt.showMM==[1,3]) + mdh = plot(ah,opt.xValues+xOffset,md,'sg','MarkerSize',12); + end + if opt.showMM == 4 + mh = plot(ah,opt.xValues+xOffset,m,'+r','Color','r','MarkerSize',12); + mdh = myErrorbar(ah,opt.xValues+xOffset,m,sem); + end + if opt.showMM == 5 + mh = plot(ah,opt.xValues+xOffset,m,'+r','Color','r','MarkerSize',12); + mdh = myErrorbar(ah,opt.xValues+xOffset,m,sd); + end + if opt.showMM == 6 + mdh(1,:) = plot(ah,squeeze(md(:,1,2:3))',repmat(md(:,1,1)',2,1),'color','r','lineWidth',2);%,'lineStyle','--'); + mdh(2,:) = plot(ah,squeeze(md(:,2,2:3))',repmat(md(:,2,1)',2,1),'color','r','lineWidth',1);%,'lineStyle','--'); + mdh(3,:) = plot(ah,squeeze(md(:,3,2:3))',repmat(md(:,3,1)',2,1),'color','r','lineWidth',1);%,'lineStyle','--'); + end + case 'flipped' + if any(opt.showMM==[1,2]) + mh = plot(ah,m,opt.xValues+xOffset,'+r','Color','r','MarkerSize',12); + end + if any(opt.showMM==[1,3]) + mdh = plot(ah,md,opt.xValues+xOffset,'sg','MarkerSize',12); + end + if opt.showMM == 4 + mh = plot(ah,m,opt.xValues+xOffset,'+r','Color','r','MarkerSize',12); + mdh = myErrorbar(ah,m,opt.xValues+xOffset,[sem,NaN(size(sem))]); + end + if opt.showMM == 5 + mh = plot(ah,m,opt.xValues+xOffset,'+r','Color','r','MarkerSize',12); + mdh = myErrorbar(ah,m,opt.xValues+xOffset,[sd,NaN(size(sd))]); + end + if opt.showMM == 6 + mdh(1,:) = plot(ah,repmat(md(:,1,1)',2,1),squeeze(md(:,1,2:3))','color','r','lineWidth',2);%,'lineStyle','--'); + mdh(2,:) = plot(ah,repmat(md(:,2,1)',2,1),squeeze(md(:,2,2:3))','color','r','lineWidth',1);%,'lineStyle','--'); + mdh(3,:) = plot(ah,repmat(md(:,3,1)',2,1),squeeze(md(:,3,2:3))','color','r','lineWidth',1);%,'lineStyle','--'); + end + end +end + +% find extents of x-axis (or y-axis, if flipped) +minX = min(opt.xValues)-stdWidth; +maxX = max(opt.xValues)+stdWidth; + +if ~isnan(xAxLim(1)) + % we have previous limits + switch opt.xyOri + case 'normal' + minX = min(minX,xAxLim(1)); + maxX = max(maxX,xAxLim(2)); + case 'flipped' + minX = min(minX,yAxLim(1)); + maxX = max(maxX,yAxLim(2)); + end +end + + +% if ~empty, use xNames +switch opt.xyOri + case 'normal' + switch opt.xMode + case 'manual' + if newAx == false + warning('DISTRIBUTIONPLOT:ERASINGLABELS','Plotting into an existing axes and specifying labels will erase previous labels') + end + set(ah,'XTick',opt.xValues); + if ~isempty(opt.xNames) + set(ah,'XTickLabel',opt.xNames) + end + case 'auto' + % no need to do anything + end + if ~isempty(opt.yLabel) + ylabel(ah,opt.yLabel); + end + % have plot start/end properly + xlim(ah,[minX,maxX]) + case 'flipped' + switch opt.xMode + case 'manual' + if newAx == false + warning('DISTRIBUTIONPLOT:ERASINGLABELS','Plotting into an existing axes and specifying labels will erase previous labels') + end + set(ah,'YTick',opt.xValues); + if ~isempty(opt.xNames) + set(ah,'YTickLabel',opt.xNames) + end + case 'auto' + % no need to do anything + end + if ~isempty(opt.yLabel) + xlabel(ah,opt.yLabel); + end + % have plot start/end properly + ylim(ah,[minX,maxX]) +end + + +%========================== + + +%========================== +%% CLEANUP & ASSIGN OUTPUT +%========================== + +if nargout > 0 + handles{1} = hh; + handles{2} = [mh;mdh]; + handles{3} = ah; + handles{4} = sh; +end + +set(ah,'NextPlot',holdState); \ No newline at end of file diff --git a/Plotting/plot_slice.m b/Plotting/plot_slice.m new file mode 100755 index 0000000..cb4f2c8 --- /dev/null +++ b/Plotting/plot_slice.m @@ -0,0 +1,88 @@ +%% This function plots the CAPs obtained for a given case +% Inputs: +% - Cp: a matrix of size n_clusters x n_voxels with the CAPs (or patterns) +% to plot +% - T: the threshold below which voxels will not be colored +% - maxC: the maximal value at which the color display will saturate +% - mask: a very long vector of logicals symbolizing the regions of the +% brain that are actually to be considered (in-brain voxels typically) +% - brain_final: a 3D matrix used to plot the greyscale brain template on which +% to overlay activity patterns +% - ai: the nii data related to the considered seed, including information +% on the scale differences between Cp indexes ('distance between voxels of +% the matrix') and actual MNI space +% - Dimension: the type of slice to plot ('X', 'Y' or 'Z') +% - MNI: the MNI coordinate of the slice to plot +% - Handle: the handle of the graph to update +function [Handle] = plot_slice(Cp,T,maxC,mask,brain_final,ai,Dimension,MNI,Handle) + + % Computes the matrix index matching the MNI coordinates of interest. + Map = inv(ai.mat); + + switch Dimension + case 'X' + ctr = round(Map(1,1)*MNI+Map(1,4)); + case 'Y' + ctr = round(Map(2,2)*MNI+Map(2,4)); + case 'Z' + ctr = round(Map(3,3)*MNI+Map(3,4)); + end + + % temp contains the volume values (within mask), and is a 3D volume after + % those lines + temp = nan(size(mask)); + temp(mask) = Cp; + temp(isnan(temp)) = 0; + temp = reshape(temp,ai.dim); + + % ho contains the structural underlay slice to plot at right slice, + % while tmpp contains the values to plot on top + switch Dimension + case {'X'} + tmpp = squeeze(temp(ctr,:,:)); + ho = squeeze(brain_final(ctr,:,:)); + case {'Y'} + tmpp = squeeze(temp(:,ctr,:)); + ho = squeeze(brain_final(:,ctr,:)); + case {'Z'} + tmpp = squeeze(temp(:,:,ctr)); + ho = squeeze(brain_final(:,:,ctr)); + end + + % I is an image with values from 0 to 1 (because the original + % image had no negative value) + I = double(ho)'/max(double(ho(:))); + I(I==0) = 1; + + % Creates an 'image with three values per pixel' + Irgb = cat(3,I,I,I); + + % Actual plotting + imagesc(Irgb,'Parent',Handle); + hold(Handle,'on'); + + % Plots the slice of interest on top of the brain template + h=imagesc(tmpp','Parent',Handle); + set(Handle,'YDir','normal'); + + % At this stage, ddd contains a set of colors, with white if the values + % are too low. We ask the axes of interest to be using this colormap + % colormap(Handle,ddd); + tmp_cm = cbrewer('div','RdBu',1000); + colormap(Handle,flipud(tmp_cm)); + + % Defines that the topmost and bottommost elements of the + % colormap will map maxC and -maxC respectively + caxis(Handle,[-1 1]*maxC); + + % Opacity: sets data points below the value of interest as + % transparent (they are white, but transparent); note that we do this + % specifically for the h imagesc (the one to plot on top) + A = ones(size(tmpp)); + A(abs(tmpp)