Update makeFaceVertexCData.m

ExampleSurfacePlotFunction.m

@@ -39,7 +39,8 @@
 end
 
 if nargin < 6 
-    climits = [nanmin(data) nanmax(data)];
+    climits = [nanmin(data(:,1)) nanmax(data(:,1))]; 
+    % use first col in case data is in annot file format
 end
 
 if nargin < 7

bluewhitered.m

@@ -0,0 +1,101 @@
+function newmap = bluewhitered(m)
+%BLUEWHITERED   Blue, white, and red color map.
+%   BLUEWHITERED(M) returns an M-by-3 matrix containing a blue to white
+%   to red colormap, with white corresponding to the CAXIS value closest
+%   to zero.  This colormap is most useful for images and surface plots
+%   with positive and negative values.  BLUEWHITERED, by itself, is the
+%   same length as the current colormap.
+%
+%   Examples:
+%   ------------------------------
+%   figure
+%   imagesc(peaks(250));
+%   colormap(bluewhitered(256)), colorbar
+%
+%   figure
+%   imagesc(peaks(250), [0 8])
+%   colormap(bluewhitered), colorbar
+%
+%   figure
+%   imagesc(peaks(250), [-6 0])
+%   colormap(bluewhitered), colorbar
+%
+%   figure
+%   surf(peaks)
+%   colormap(bluewhitered)
+%   axis tight
+%
+%   See also HSV, HOT, COOL, BONE, COPPER, PINK, FLAG, 
+%   COLORMAP, RGBPLOT.
+if nargin < 1
+   m = size(get(gcf,'colormap'),1);
+end
+bottom = [0 0 0.5];
+botmiddle = [0 0.5 1];
+middle = [1 1 1];
+topmiddle = [1 0 0];
+top = [0.5 0 0];
+% Find middle
+lims = get(gca, 'CLim');
+% Find ratio of negative to positive
+if (lims(1) < 0) & (lims(2) > 0)
+    % It has both negative and positive
+    % Find ratio of negative to positive
+    ratio = abs(lims(1)) / (abs(lims(1)) + lims(2));
+    neglen = round(m*ratio);
+    poslen = m - neglen;
+
+    % Just negative
+    new = [bottom; botmiddle; middle];
+    len = length(new);
+    oldsteps = linspace(0, 1, len);
+    newsteps = linspace(0, 1, neglen);
+    newmap1 = zeros(neglen, 3);
+
+    for i=1:3
+        % Interpolate over RGB spaces of colormap
+        newmap1(:,i) = min(max(interp1(oldsteps, new(:,i), newsteps)', 0), 1);
+    end
+
+    % Just positive
+    new = [middle; topmiddle; top];
+    len = length(new);
+    oldsteps = linspace(0, 1, len);
+    newsteps = linspace(0, 1, poslen);
+    newmap = zeros(poslen, 3);
+
+    for i=1:3
+        % Interpolate over RGB spaces of colormap
+        newmap(:,i) = min(max(interp1(oldsteps, new(:,i), newsteps)', 0), 1);
+    end
+
+    % And put 'em together
+    newmap = [newmap1; newmap];
+
+elseif lims(1) >= 0
+    % Just positive
+    new = [middle; topmiddle; top];
+    len = length(new);
+    oldsteps = linspace(0, 1, len);
+    newsteps = linspace(0, 1, m);
+    newmap = zeros(m, 3);
+
+    for i=1:3
+        % Interpolate over RGB spaces of colormap
+        newmap(:,i) = min(max(interp1(oldsteps, new(:,i), newsteps)', 0), 1);
+    end
+
+else
+    % Just negative
+    new = [bottom; botmiddle; middle];
+    len = length(new);
+    oldsteps = linspace(0, 1, len);
+    newsteps = linspace(0, 1, m);
+    newmap = zeros(m, 3);
+
+    for i=1:3
+        % Interpolate over RGB spaces of colormap
+        newmap(:,i) = min(max(interp1(oldsteps, new(:,i), newsteps)', 0), 1);
+    end
+
+end

checkVertsFacesRoisData.m

@@ -0,0 +1,133 @@
+function [verts,faces,rois,data] = checkVertsFacesRoisData(varargin)
+%% Checks size and shape of data describing brain surface mesh
+% Ensures that verts is V x 3, faces is F x 3, rois is V x 1, and data is V x 1
+% or R x 1
+%
+%
+%% Syntax
+%  [verts,faces,rois,data] = checkVertsFacesRoisData(verts,faces,rois,data)
+%
+%
+%% Input Arguments
+%  verts - xyz coordinates of vertices (three column matrix)
+%  faces - triangulation ie. IDs of vertices making up each face (three column matrix) 
+%  rois - roi allocation of each vertex (V x 1 vector)
+%  data - data allocated to each vertex or roi (V x 1 or R x 1 vector)
+%
+%
+%% Name-Value Arguments
+%  checkContents - flag to check values within matrices (true (default) | false)
+%  If `checkContents` is set to false, only the shape of the input matrices will
+%  be tested. If `checkContents` is true, then the contents will also be checked
+%  e.g. that `faces` uses all the vertices and that `rois` cannot have negative
+%  values. 
+%
+%  fillEmpty - flag to populate `rois` and `data`, if the inputs are empty (false (default) | true) 
+%  If set to true, rois will be set to all ones, and data will be set to the
+%  second column of `verts`.
+%
+%
+%% Output Arguments
+%  verts, faces, rois, data - transposed if needed
+% 
+% 
+%% See Also
+% plotSurfaceROIBoundary, read_vtk
+% 
+% 
+%% Authors 
+% Mehul Gajwani, Monash University, 2023
+% 
+% 
+
+
+%% Prelims
+ip = inputParser;
+validationFcn = @(x) (isnumeric(x) && ismatrix(x));
+addOptional(ip, 'verts', [], @(x) validationFcn(x));
+addOptional(ip, 'faces', [], @(x) validationFcn(x));
+addOptional(ip, 'rois',  [], @(x) validationFcn(x) || islogical(x));
+addOptional(ip, 'data',  [], @(x) validationFcn(x) || islogical(x));
+
+addParameter(ip, 'checkContents', true, @islogical);
+addParameter(ip, 'fillEmpty', false, @islogical);
+
+parse(ip, varargin{:});
+verts = ip.Results.verts;
+faces = ip.Results.faces;
+rois  = +ip.Results.rois;
+data  = +ip.Results.data;
+
+
+%% Check shape
+% Check second dimension and transpose if needed
+
+matrices = {verts, faces, rois, data};
+names = {'Vertices', 'Faces', 'ROIs', 'Data'};
+target2ndDim = {3, 3, 1, 1};
+
+for ii = 1:length(matrices)
+    if isempty(matrices{ii})
+        % skip - do nothing
+    elseif size(matrices{ii}, 2) == target2ndDim{ii}
+        % good - do nothing
+    elseif size(matrices{ii}, 1) == target2ndDim{ii}
+        matrices{ii} = matrices{ii}.';
+    else
+        error('%s must have a dimension of size %d', names{ii}, target2ndDim{ii});
+    end
+end
+
+[verts, faces, rois, data] = deal(matrices{:});
+
+if isallhere({verts, rois})
+    assert(size(verts, 1) == size(rois, 1), ...
+        'Vertices and ROIs must have the same number of points');
+end
+
+
+%% Fill empty, if desired
+if ip.Results.fillEmpty
+    if isempty(rois); rois = ones(size(verts, 1),1); end
+    if isempty(data); data = verts(:, 2); end
+end
+
+
+%% Check contents
+if ip.Results.checkContents
+    if isallhere({verts, faces})
+        if max(faces, [], "all") ~= size(verts, 1)
+            warning('The vertices are not all mapped to the faces');
+        end
+    end
+
+    if isallhere({faces}); mustBeNonnegative(faces); mustBeInteger(faces); end
+
+    if isallhere({rois})
+        mustBeNonnegative(rois); mustBeInteger(rois);
+
+        goodRois = ( all(rois) && (length(unique(rois))==max(rois)) ) || ... % no non-zeros rois
+            (~all(rois) && (length(unique(rois))==max(rois)+1) ); % some non-zero rois
+
+        if ~goodRois
+            warning('ROIs appear to not be sequential');
+        end
+    end
+
+    if isallhere({data, verts}) || isallhere({data, rois})
+        assert(size(data, 1) == size(verts, 1) || size(data, 1) == max(rois), ...
+            'Data should be one per vertex or one per ROI');
+    end
+end
+
+end % main
+
+function out = isallhere(inp)
+% returns true if inp is (i) a non-empty matrix OR (ii) a cell with ALL non-empty entries
+% else returns false
+if ~iscell(inp); out = ~isempty(inp); return; end
+out = ~any(cellfun(@isempty, inp));
+end
+
+
+

demo_plotSurfaceROIBoundary.m

@@ -16,7 +16,7 @@
 plotSurfaceROIBoundary(surface,lh_rand200,1:100,'faces',jet(100),2);
 
 % The following options set up the patch object to look pretty. This works
-% well for the left hemisphere (medial and lateral). Change the inputs to 
+% well for the left hemisphere (medial and lateral). Change the inputs to
 % 'view' to see the brain from different angles ([-90 0] for left and [90 0]
 % for right I find works well)
 
@@ -26,8 +26,7 @@
 axis off
 axis tight
 axis equal
-axis vis3d
-
+%%
 ax2 = axes('Position',[0.01+(1/3) 0 .3 1]);
 
 cmap = flipud(hot(130));
@@ -45,14 +44,13 @@
 axis off
 axis tight
 axis equal
-axis vis3d
-
+%%
 ax3 = axes('Position',[0.01+(2/3) 0 .3 1]);
 
 % This plots sulcal depth, which is defined for each vertex
 
 surface.vertices = lh_verts;
-plotSurfaceROIBoundary(surface,lh_aparc,lh_sulc,'centroid',parula(100),4);
+plotSurfaceROIBoundary(surface,lh_aparc,lh_sulc,'centroid',parula(100),1);
 camlight(80,-10);
 camlight(-80,-10);
 
@@ -63,15 +61,15 @@
 axis equal
 axis vis3d
 
-% Demonstrate different types of plots
+%% Demonstrate different types of plots
 
 surface.vertices = lh_inflated_verts;
 boundary_type = {'faces','midpoint','centroid','edge_vertices','edge_faces',...
     'faces','midpoint','centroid','edge_vertices','edge_faces'};
 linewidth = 4;
 
 
-% Set up the colors for each vertex. This mirrors how the redone colour 
+% Set up the colors for each vertex. This mirrors how the redone colour
 % assignment is performed by makeFaceVertexCData
  lh_rand200_color_map = lines(34);
  lh_rand200_ = lh_rand200;
@@ -82,11 +80,11 @@
  lh_rand200_color(isnan(lh_rand200_),:) = .5;
 
  FaceVertexCData = makeFaceVertexCData(surface.vertices,surface.faces,lh_rand200,lh_rand200,lh_rand200_color_map);
- 
+
 for i = 1:10
 
 figure('Position',[0 0  1680 933])
-    
+
     % The data here is just each ROIs own ID number
 
         if i < 6
@@ -99,10 +97,10 @@
             if i == 6
             savename = ['sulc_',boundary_type{i},'_flat.png'];
             else
-            savename = ['sulc_',boundary_type{i},'_interp.png'];    
+            savename = ['sulc_',boundary_type{i},'_interp.png'];
             end
         end
-    
+
     p = plotSurfaceROIBoundary(surface,lh_rand200,data,boundary_type{i},cmap,linewidth);
 
     camlight(80,-10);
@@ -114,28 +112,28 @@
     axis tight
     axis equal
     %axis vis3d
-    
+
     % Mapping on the ROI id of each vertex to help with understanding how
     % this all works
     hold on
-        
+
     s = scatter3(lh_inflated_verts(:,1)*1.01,lh_inflated_verts(:,2),lh_inflated_verts(:,3),40,lh_rand200_color,'filled');
     s.Clipping = 'off';
     s.MarkerEdgeColor = 'k';
     s.LineWidth = 1;
-    
-    % This just zooms into the area of interest    
+
+    % This just zooms into the area of interest
     ylim([-25.2699   -8.7600])
     zlim([20.2174   31.3705])
 
     p.EdgeColor = 'k';
     p.EdgeAlpha = .5;
-    
+
     % If you wanted to make a colorbar, this is what you would have to do:
     % colormap(cmap)
     % caxis([min(data) max(data)])
     % c = colorbar
-    
+
     %print(['./figures/',savename],'-dpng')
 
 end