[Master Index]
[Index for Toolbox]
bst_imaging_display
(Toolbox/bst_imaging_display.m in BrainStorm 2.0 (Alpha))
Function Synopsis
varargout = bst_imaging_display(varargin)
Help Text
BST_IMAGING_DISPLAY - : display utility for cortical distributed source imaging
function varargout = bst_imaging_display(varargin)
function varargout = bst_imaging_display(ImagingFigure,ImageGridLoc,ImageGridAmp,varargin)
if nargin == 2
ImageGridAmp = varargin{1};
OPTIONS = varargin{2};
----------------------------------------------------------------------------------------------
Cross-Reference Information
This function calls
This function is called by
Listing of function C:\BrainStorm_2001\Toolbox\bst_imaging_display.m
function varargout = bst_imaging_display(varargin)
%BST_IMAGING_DISPLAY - : display utility for cortical distributed source imaging
% function varargout = bst_imaging_display(varargin)
%function varargout = bst_imaging_display(ImagingFigure,ImageGridLoc,ImageGridAmp,varargin)
% if nargin == 2
% ImageGridAmp = varargin{1};
% OPTIONS = varargin{2};
%
% ----------------------------------------------------------------------------------------------
%<autobegin> ---------------------- 12-Oct-2004 01:10:59 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Visualization
%
% Alphabetical list of external functions (non-Matlab):
% toolbox\bst_message_window.m
% toolbox\catci.m
% toolbox\curvature_cortex.m
% toolbox\dataplot_cb.m
% toolbox\get_user_directory.m
% toolbox\good_channel.m
% toolbox\grayish.m
% toolbox\interp_mail.m
% toolbox\norcol.m
% toolbox\norlig.m
% toolbox\plot_dipole.m
% toolbox\vertices_connectivity.m
% toolbox\view_surface.m
%
% Subfunctions in this file, in order of occurrence in file:
% [] = DataOnSensors % Sensors are color-coded according to the data at each location
% [] = DataOnScalp % Interpolated Scalp Topography of Surface Data
% [] = CorticalMap(TessWin,OPTIONS) % See cortical current maps interpolated on the proper cortical surface
% [mixedRGB] = mix_curvature_current(varargin)
%
% Figure Files opened by this function:
% 'mapping.fig'
% 'tessellation_window.fig'
%
% Format of strings below: Type:Style:Tag, "String", CallBack Type and Call
% <automatic> callback is <Tag>_Callback by Matlab default
%
% Callbacks by figure mapping.fig
% figure::mapping "" uses ButtonDownFcn for dataplot_cb line_select
% uicontrol:checkbox:ABSOLUTE "Absolute" uses Callback for dataplot_cb mutincomp_cmapping
% uicontrol:checkbox:Colorbar "Colorbar" uses Callback for dataplot_cb Colorbar_mapping
% uicontrol:checkbox:fit "Sensor array" uses Callback for dataplot_cb mapping_head_shape
% uicontrol:checkbox:gradient "Gradient" uses Callback for dataplot_cb mapping_type
% uicontrol:checkbox:magnetic "Mag. Field Estim." uses Callback for dataplot_cb mapping_type
% uicontrol:checkbox:movie "Movie" uses Callback for dataplot_cb mapping_display_type
% uicontrol:checkbox:raw "Raw Mapping" uses Callback for dataplot_cb mapping_type
% uicontrol:checkbox:RELATIVE "Relative" uses Callback for dataplot_cb mutincomp_cmapping
% uicontrol:checkbox:scalp "Scalp surface" uses Callback for dataplot_cb mapping_head_shape
% uicontrol:checkbox:ShowSensorLabels "Sensor Labels" uses Callback for dataplot_cb mapping_type
% uicontrol:checkbox:single "Single" uses Callback for dataplot_cb mapping_display_type
% uicontrol:checkbox:slides "Slides" uses Callback for dataplot_cb mapping_display_type
% uicontrol:checkbox:sphere "Sphere" uses Callback for dataplot_cb mapping_head_shape
% uicontrol:edit:current_time "0" uses Callback for dataplot_cb set_current_time
% uicontrol:pushbutton:go "Apply" uses Callback for dataplot_cb gomapping
% uicontrol:pushbutton:Pushbutton1 "Quit" uses Callback for delete(gcbf)
% uicontrol:slider:slider_time "" uses Callback for dataplot_cb mapping_slider
%
% Callbacks by figure tessellation_window.fig
% figure::tessellation_window "" uses CreateFcn for movegui northwest
%
% At Check-in: $Author: Mosher $ $Revision: 28 $ $Date: 10/11/04 11:32p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 23-Sep-2004
%
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
% University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
% Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
% CNRS, Hopital de la Salpetriere, Paris, France
%
% See BrainStorm website at http://neuroimage.usc.edu for further information.
%
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 12-Oct-2004 01:10:59 -----------------------
% EK 4-Oct-2004 Minor change in displaying the cortex: If the
% ImageGridAmp is a vector (like in Neural Index Case) then
% display this for any time point selected in time series window.
% ----------------------------------------------------------------------------
if nargin == 4
OPTIONS = varargin{1};
end
Users = get_user_directory;
% Define default values for display options
Def_CutPlane_OPTIONS = struct(...
'Top', 0, ...
'Bottom', 0, ...
'Front', 0, ...
'Back', 0,...
'Right', 0,...
'Left', 0,...
'X', [],...
'Y', [],...
'Z', [] ...
);
Def_OPTIONS = struct(...
'CurvatureColor1',[.4 .4 .4],...
'CurvatureColor2',[.6 .6 .6],...
'Dipoles',[],...
'OrthoViews',0,...
'Time', 1, ...
'ScalpMaps',0,...
'ColorSensors',0,...
'CorticalMap',1,...
'ColorMapping','normalized',...
'CDThreshold',.5,...
'currDensityTransparency',.2, ...
'CurvatureThreshold',.16,...
'FaceVertexAlpha',.2,...
'hPatch', NaN,...
'MapCurvature',1,...
'timeText',NaN,...
'timeValue',[],...
'Curvature', [], ...
'InitialCurvature',[],...
'FaceAlpha',1,...
'VertConn', [],...
'CutPlane', Def_CutPlane_OPTIONS...
);
% Check arguments passed to and required from function call
if nargin < 3 % return default OPTIONS values
if nargout > 1
varargout{1} = Def_OPTIONS;
varargout{2} = Def_OPTIONS;
else
varargout{1} = Def_OPTIONS;
end
if nargin == 2
ImageGridAmp = varargin{1};
if isempty(ImageGridAmp)
ImageGridAmp = 0; % no current mapping - just a look at anatomy
end
OPTIONS = varargin{2};
else
return
end
else
ImagingFigure = varargin{1};
ImageGrid = varargin{2};
ImageGridAmp = varargin{3};
OPTIONS = varargin{4};
end
% Check field names of passed OPTIONS and fill missing ones with default values
DefFieldNames = fieldnames(Def_OPTIONS);
for k = 1:length(DefFieldNames)
if ~isfield(OPTIONS,DefFieldNames{k})
OPTIONS = setfield(OPTIONS,DefFieldNames{k},getfield(Def_OPTIONS,DefFieldNames{k}));
end
end
clear Def_OPTIONS
% ----------------------------------------------------------------------------------------------
% tmp = findobj(0,'tag','subplot_tess_window');
% if ~isempty(tmp)
% ImagingFigure = tmp; % Consider the slides or movie window
% end
if ~ishandle(OPTIONS.hPatch) % no surface available: load vertex/faces information
% Load surface information ------------------
try
load(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1})) % CBB - suppose there is one surface per tessellation file
catch
try
load(ImageGrid.GridLoc{1})
catch
errordlg(sprintf('Subject file %s could not be found. Please update subject''s entry in database or manually update the subject information in selected results file.',fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1})),'Unvalid subject database entry');
return
end
end
FigName = Comment{ImageGrid.iGrid(1)}; clear Comment % Name of figure where to plot the surface
%OPTIONS.SurfaceName = FigName;
fv.faces = Faces{ImageGrid.iGrid(1)}; clear Faces
fv.vertices = Vertices{ImageGrid.iGrid(1)}'; clear Vertices
nverts_max = 10000; % Max number of vertices before warning
if size(fv.vertices,1) > 10000
ButtonName=questdlg(sprintf('Surface has %d vertices; would you like to downsample this number ?',size(fv.vertices,1)), ...
'Large Surface Detected', ...
'Yes','No','Yes');
switch(ButtonName)
case 'Yes'
hmsg = msgbox('Downsampling original surface. . .','Working!','warn'); drawnow
fv = reducepatch(fv,10000/size(fv.vertices,1),'fast','verbose');
delete(hmsg);
otherwise
% do nothing
end
end
% Load surface information----DONE----------
figure(ImagingFigure)
% Figure Parameters ------------------------
set(ImagingFigure,'CurrentAxes',findobj(ImagingFigure,'Tag','MainAxis'),...
'Renderer','OpenGL',...
'Color','w', 'Colormap',gray,'Menubar','figure','Name',FigName);
else
fv.vertices = get(OPTIONS.hPatch,'Vertices');
end
% % Surface orientation
% if ~isfield(Visu.Data,'System')
% Visu.Data.System = 'ctf';
% end
% Figure Parameters -----------DONE---------
if 0 % deprecated code (Sb 08-Oct-2003)
%FragmentMenu = findobj(gcbf,'Tag','FragmentMenu');
%nclust = get(FragmentMenu,'Value')-1; % Get label of clustering to visualize
nclust = 0;
hpatch = []; % handle(s) to patch objects in figure
if OPTIONS.OrthoViews == 0
hpatch = findobj(findobj(ImagingFigure,'Tag','MainAxes'),'Type','patch','Visible','on'); % CBB - store these handles in figure appdata ? (SB)
else
if ~isempty(findobj(findobj(ImagingFigure,'Tag','axsub1'),'Type','patch'))
hpatch(1) = findobj(findobj(ImagingFigure,'Tag','axsub1'),'type','patch');
hpatch(2) = findobj(findobj(ImagingFigure,'Tag','axsub2'),'type','patch');
hpatch(3) = findobj(findobj(ImagingFigure,'Tag','axsub3'),'type','patch');
hpatch(4) = findobj(findobj(ImagingFigure,'Tag','axsub4'),'type','patch');
else
if isempty(hpatch)
hpatch = [];
end
end
end
end % Deprecated code -------------------------
% Work on cut planes --------------------------
% work with facelapha property only?
% Work on cut planes -------------DONE---------
% Compute or retrieve curvature information --
if OPTIONS.MapCurvature % Color coding of surface curvature
if isempty(OPTIONS.Curvature) % No curvature mapping is available - compute it
% Is vertex connectivity available ?
%load(fullfile(Users.SUBJECTS,Users.CurrentData.SubjectFile),'VertConn');
try
load(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1}),'VertConn');
catch
ImageGrid.GridLoc{1} = strrep(ImageGrid.GridLoc{1},Users.SUBJECTS,'');
end
load(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1}),'VertConn');
if ~exist('VertConn','var') % Vertex connectivity not computed yet
load(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1}),'Comment');
VertConn = cell(length(Comment),1);
end
%VertConnFile = VertConn; % Filename containing connectivity information
[pathname, filename] = fileparts(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1}));
cd(pathname)
%if ~exist(VertConnFile,'file')
if isempty(VertConn{ImageGrid.iGrid(1)})
% How many surfaces in tessellation file ?
% should be only one but older BrainStorm files may contain several
load(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1}),'Comment');
nsurfaces = length(Comment); % Number of surfaces in tessellation file
VertConn = cell(nsurfaces,1); % Create cell array where to store vertex connectivity values
VertConn{ImageGrid.iGrid(1)} = vertices_connectivity(fv); % compute vertex connectivity
[pathname, filename] = fileparts(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1}));
cd(pathname)
save(ImageGrid.GridLoc{1},'VertConn','-append');
else % load existing vertex connectivity map
%load(fullfile(pathname,VertConnFile))
end
% Now compute curvature information
OPTIONS.VertConn = VertConn{ImageGrid.iGrid(1)}; clear VertConn
%[vertexcolorh,map]=c_activ(FV,TessWinUserData.VertConn); clear FV ; OPTIONS.Curvature = vertexcolorh;
% Parameters for curvature computation
sigmoid_const = .1;
show_sigmoid = 0;
if isempty(OPTIONS.InitialCurvature)
load(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1}),'Curvature');
if ~exist('Curvature','var') % Vertex connectivity not computed yet
load(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1}),'Comment');
Curvature = cell(length(Comment),1);
end
if isempty(Curvature{ImageGrid.iGrid(1)}) | size(Curvature{ImageGrid.iGrid(1)},2) > 1 % Older curvature computation
bst_message_window('wrap',...
{' ','Computing surface curvature. . .'})
[curvature_sigmoid,Curvature{ImageGrid.iGrid(1)}] = curvature_cortex(fv,OPTIONS.VertConn,sigmoid_const,show_sigmoid);
bst_message_window('wrap',...
{'Computing surface curvature ->DONE',' '})
clear curvature_sigmoid
save(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1}),'Curvature','-append');
end
end
OPTIONS.InitialCurvature = Curvature{ImageGrid.iGrid(1)}; clear Curvature
end
% apply threshold
OPTIONS.Curvature(OPTIONS.InitialCurvature >= OPTIONS.CurvatureThreshold) = .8;
OPTIONS.Curvature(OPTIONS.InitialCurvature < OPTIONS.CurvatureThreshold) = -.8;
end
% Compute or retrieve curvature information --
% Colormap for current density---
tmp = hot(140);%flipdim(hot(100),1);
currDensityColorMap = tmp(1:100,:); % Current density colormap
%anatColorMap = repmat(linspace(.4,.6,100)',1,3); % Anatomy colormap
anatColorMap = [OPTIONS.CurvatureColor1;...
zeros(98,3);...
OPTIONS.CurvatureColor2]; % Anatomy colormap - articifical zero padding
% Colormap ----------DONE--------
if ImageGridAmp == 0 % No activation is specified - look at anatomy only
ImageGridAmp = ones(size(fv.vertices,1),length(OPTIONS.Time));
end
% Create patch object -------------------------
if ~ishandle(OPTIONS.hPatch) % No patch handle is passed (i.e. need to create one)
figure(ImagingFigure)
PlotHandles.SurfaceName = FigName;
[PlotHandles.hFigure,PlotHandles.hPatch,PlotHandles.hLights] = ...
view_surface(FigName,fv.faces,fv.vertices,ImageGridAmp(:,OPTIONS.Time));
% Dull rendering
delete(findobj(PlotHandles.hFigure,'type','light'))
%view(-10,0),
material dull
%PlotHandles.hLights = camlight;
set(PlotHandles.hPatch,'backfacelighting','unlit')
set(gca,'AmbientLightColor',[1 1 1])
hold on
rotate3d(PlotHandles.hFigure,'on')
%set(hf,'color','k')
else
PlotHandles.hPatch = OPTIONS.hPatch;
end
% Create patch object ----------DONE-----------
% Translate anatomy curvature in RGB values
if ~OPTIONS.MapCurvature % no color coding of surface curvature
if ishandle(PlotHandles.hPatch)
fv.vertices = get(PlotHandles.hPatch,'Vertices');
end
index_anatomy = 100*ones(size(fv.vertices,1),1);
else
min_curv = min(OPTIONS.Curvature);
max_curv = max(OPTIONS.Curvature);
index_anatomy = floor((OPTIONS.Curvature-min_curv)*99/(max_curv-min_curv)+1);
end
anatRGB = anatColorMap(index_anatomy,:);
% Translate current density values in RGB
% Amplitude scaling
ImageGridAmp = double(ImageGridAmp);
switch lower(OPTIONS.ColorMapping) % is color mapping time-relative or time-normalized ?
case 'normalized'
if size(ImageGridAmp,2) == 1 %E.g Representation of neural activity index
if OPTIONS.Time >1
bst_message_window('wrap',{'Currently visualized result does not change with time',' '});
end
currDensity = abs(ImageGridAmp(:,1))/max(abs(ImageGridAmp(:,1)));
else
currDensity = abs(ImageGridAmp(:,OPTIONS.Time))/max(abs(ImageGridAmp(:,OPTIONS.Time)));
end
limitValue = max(currDensity); % Store maximum value
case 'relative'
if size(ImageGridAmp,2) == 1 %E.g Representation of neural activity index
if OPTIONS.Time >1
bst_message_window('wrap',{'Currently visualized result does not change with time',' '});
end
currDensity = abs(ImageGridAmp(:,1))/max(abs(ImageGridAmp(:)));
else
currDensity = abs(ImageGridAmp(:,OPTIONS.Time))/max(abs(ImageGridAmp(:)));
end
limitValue = 1; % Force it to 1
end
% Apply amplitude threshold
%tmpcurrDensity = currDensity;
currDensity(currDensity < OPTIONS.CDThreshold) = 0;
if limitValue == min(currDensity) % no current value above threshold
limitValue = -1 ;
end
%clear tmpcurrDensity
if 1%(viewmode == 1)%if real values
index_density = floor((currDensity + limitValue)*99/(2*limitValue)+1);
else %if abs values
currDensity = abs(currDensity);
index_density = floor((currDensity-limitValue)*99/limitValue+100);
end
currDensityRGB = currDensityColorMap(index_density,:);
% Now mix curaveture and current RGBs
mixedRGB = anatRGB;
toBlend = find(currDensity ~= 0); % Find vertex indices holding non-zero activation (after thresholding)
mixedRGB(toBlend,:) = OPTIONS.currDensityTransparency*anatRGB(toBlend,:)+...
(1-OPTIONS.currDensityTransparency)*currDensityRGB(toBlend,:);
% Assign color to surface elements
set(PlotHandles.hPatch,'FaceVertexCData',mixedRGB);
drawnow
% Add text label to indicate current time
if ~isempty(OPTIONS.timeValue)
if ~ishandle(OPTIONS.timeText)
OPTIONS.timeText = text(0.05,0.05,0,sprintf('%3.1f ms',1000*OPTIONS.timeValue),'units','normalized');
set(OPTIONS.timeText,'fontweight','normal','color',[.3 .3 .3],'Fontname','helvetica','Fontsize',8,'FontUnits','Point');
else
set(OPTIONS.timeText,'String',sprintf('%3.1f ms',1000*OPTIONS.timeValue))
end
end
% Now overlay dipole models if any -------------------------------------
%rotate3d on
if isempty(OPTIONS.Dipoles)
varargout{1} = PlotHandles;
varargout{2} = OPTIONS;
%rotate3d(get(PlotHandles.hPatch,'Parent'),'on')
return
end
%set(PlotHandles.hPatch,'edgecolor','none','FaceAlpha',0.3); hold on
% Make transparency decrease with depth
%fv.vertices = fv.vertices - repmat( mean(fv.vertices),size(fv.vertices,1),1);
% [TH,PHI,R] = CART2SPH(fv.vertices(:,1),fv.vertices(:,2),fv.vertices(:,3)); clear TH PHI
% R = R/max(R);
% set(PlotHandles.hPatch,'edgecolor','none','FaceVertexAlpha',1-R,'facealpha','flat'); hold on
axes(get(PlotHandles.hPatch,'Parent'))
set(PlotHandles.hPatch,'FaceVertexAlpha',OPTIONS.FaceVertexAlpha,'facealpha','flat'); hold on
scale = norcol(OPTIONS.Dipoles.TimeSeries); % Scale factor for amplitude display
scale = scale/max(scale);
scale(scale <.5) = .5; % Avoid too small dipoles
scale = scale/200; % Dipole size needs to fit within cortical volume
for src = 1:length(OPTIONS.Dipoles.SourceLoc)
axes(get(PlotHandles.hPatch,'Parent'))
if isempty(OPTIONS.Dipoles.Color{src})
OPTIONS.Dipoles.Color{src} = rand(1,3);
end
plot_dipole(OPTIONS.Dipoles.SourceLoc{src}',...
OPTIONS.Dipoles.Orientation{src}',...
scale(src),OPTIONS.Dipoles.Color{src});
end
varargout{1} = PlotHandles;
varargout{2} = OPTIONS;
rotate3d(PlotHandles.hFigure,'on')
return
% Deprecated code below
% Call subfunctions depending on type fo display
if OPTIONS.ColorSensors
DataOnSensors
end
if OPTIONS.ScalpMaps
DataOnScalp
end
if OPTIONS.CorticalMap
%CorticalMap(ImagingFigure,hpatch)
else
if get(OPTIONS.MapCurvature,'Value')
set(hpatch,'FaceVertexCData',vertexcolorh)
end
end
set(hpatch,'Visible','on')
vertexcolorh = get(hpatch,'FaceVertexCdata');
% if nclust > 0%k == nsurf & nclust > 0 % Some clustering has been selected - display clusters on the surface
% set(h,'FacevertexCdata',Clusters{nsurf}.Cluster{nclust},'Userdata',Clusters{nsurf}.Cluster{nclust}')
% colormap(rand(Clusters{nsurf}.NClusters(nclust),3))
% end
hold on
set(hpatch,'edgecolor','none','facealpha',OPTIONS.FaceAlpha)
% Update envelope properties ------------------------------------------------------
ctet = get(hpatch,'vertices');
siz = size(vertexcolorh,2);
if size(vertexcolorh,1) == 1
vertexcolorh = ones(size(ctet,1),1)*vertexcolorh;
end
if 0
if Vals(k,3) == 1
vertexcolorh(find(ctet(:,2) < 0),:) = NaN * ones(length(find(ctet(:,2) < 0)),siz);
end
if Vals(k,4) == 1
vertexcolorh(find(ctet(:,2) > 0),:) = NaN * ones(length(find(ctet(:,2) > 0)),siz);
end
if Vals(k,3) == 1 | Vals(k,4) == 1 % Separate both hemispheres
% centering of the brain
ctet2(:,1)=ctet(:,1)- mean(ctet(:,1));%sum(sum(ctet(:,1)))/size(ctet,1).*ones(size(ctet,1),1);
ctet2(:,2)=ctet(:,2)- mean(ctet(:,2));%sum(sum(vu(:,2)))/size(vu,1).*ones(size(vu,1),1);
ctet2(:,3)=ctet(:,3)- mean(ctet(:,3));;%sum(sum(vu(:,3)))/size(vu,1).*ones(size(vu,1),1);
% ctet2(:,1)=ctet(:,1)- sum(sum(ctet(:,1)))/size(ctet,1).*ones(size(ctet,1),1);
% ctet2(:,2)=ctet(:,2)- sum(sum(ctet(:,2)))/size(ctet,1).*ones(size(ctet,1),1);
% ctet2(:,3)=ctet(:,3)- sum(sum(ctet(:,3)))/size(ctet,1).*ones(size(ctet,1),1);
%
% % find short and long axes
% larg1=max(max(ctet2(:,1)))-min(min(ctet2(:,1)));
% larg2=max(max(ctet2(:,2)))-min(min(ctet2(:,2)));
% if larg1<larg2
% short_axe=1;
% else
% short_axe=2;
% end
%
% cd(Users.SUBJECTS)
% load(Users.CurrentData.SubjectFile,'VertConn'); % Is the vertex connectivity available ?
% load(VertConn)
% xVertConn = VertConn;
% if isempty(xVertConn{nsurf}) % Vertex connectivity not available: compute and save.
% FV.faces = get(h,'faces');
% FV.vertices = get(h,'vertices');
% tmp = vertices_connectivity(FV); clear FV
% xVertConn{nsurf} = tmp; clear tmp
% load(Users.CurrentData.SubjectFile,'VertConn');
% VertConnFile= VertConn; VertConn = xVertConn; clear xVertConn
% save(VertConnFile,'VertConn','-append')
% VertConn = VertConn{nsurf};
% else
% VertConn = xVertConn{nsurf};
% end
%
% [I_pos,I_neg] = separation(ctet2,VertConn); clear ctet2
%SVD the vertex coordinates to find principal axes
%Dowsize the vertex location matrix to run the SVD and find the axis of principal direction for the cortical volume
nverts = size(ctet2,1);
FV.faces = get(h,'faces');
FV.vertices = get(h,'vertices');
nfv = reducepatch(FV,1000); clear FV
tmp = nfv.vertices; clear nfv
tmp = tmp - repmat(mean(tmp),size(tmp,1),1);
ndown = 500; % Objective for the number of vertices considered for determination of the main axes.
rate = ceil(nverts/ndown);
tmp = ctet2(1:rate:end,:);
[U,s,V] = svd(tmp'/sqrt(size(tmp,1)-1),0); % Compute the SVD of the covariance coordinate matrix.
s = sqrt(diag(s)); % Get the standard deviation in each direction
% Draw orthogonal planes to illustrate the 3 orthogonal axes
% Brain centroid
brain_c = mean(ctet);
% Get the axes limits;
axh = get(h,'Parent'); % Get current surface parent
Xlims = get(axh,'Xlim'); Ylims = get(axh,'Ylim'); Zlims = get(axh,'Zlim');
% Parametric equation for all lines passing by brain_c in 3D: X = brain_c + U(:,i)*t
fline = inline('brain_c'' + u*x','x','brain_c','u');% Create a function for the parametric expression of the line
% Sagittal Line
p1 = fline(-.1,brain_c,U(:,1))';
p2 = fline(.1,brain_c,U(:,1))';
hline(1) = line([p1(1) p2(1)],[p1(2) p2(2)],[p1(3) p2(3)]);
% Sagittal Line
p1 = fline(-.1,brain_c,U(:,2))';
p2 = fline(.1,brain_c,U(:,2))';
hline(2) = line([p1(1) p2(1)],[p1(2) p2(2)],[p1(3) p2(3)]);
% Sagittal Line
p1 = fline(-.1,brain_c,U(:,3))';
p2 = fline(.1,brain_c,U(:,3))';
hline(3) = line([p1(1) p2(1)],[p1(2) p2(2)],[p1(3) p2(3)]);
% % Select the clustering seed in one of the hemisphere: pick up the extreme point left or right vertex
% if short_axe == 1
% [mm, imax] = min(abs(ctet2(:,1)));
% else
% [mm, imax] = min(abs(ctet2(:,2)));
% end
%
% [CLASS, NumClass] = crtx_cluster(VertConn,2,imax,1);
%
% Get corresponding t values
if Vals(k,3) == 1
I_pos = find(CLASS{1} == 1);
vertexcolorh(I_pos,:) = NaN * ones(length(I_pos),siz);
elseif Vals(k,4) == 1
I_neg = find(CLASS{1} == 2);
vertexcolorh(I_neg,:) = NaN * ones(length(I_neg),siz);
end
if Vals(k,5) == 1
vertexcolorh(find(ctet(:,1) > 0),:) = NaN * ones(length(find(ctet(:,1) > 0)),siz);
end
if Vals(k,6) == 1
vertexcolorh(find(ctet(:,1) < 0),:) = NaN * ones(length(find(ctet(:,1) < 0)),siz);
end
if Vals(k,7) == 1
vertexcolorh(find(ctet(:,3) > 0),:) = NaN * ones(length(find(ctet(:,3) > 0)),siz);
end
if Vals(k,8) == 1
vertexcolorh(find(ctet(:,3) < 0),:) = NaN * ones(length(find(ctet(:,3) < 0)),siz);
end
end
end
clear ctet
if iscell(vertexcolorh)
set(h,'FaceVertexCData',vertexcolorh{1},'Visible','on')
if length(hpatch)>1 % OrthoViews
for kk = 1:length(hpatch)
h = hpatch(kk);
set(h,'FaceVertexCData',vertexcolorh{1},'Visible','on')
end
end
else
if OPTIONS.OrthoViews == 1
set(hpatch,'FaceVertexCData',vertexcolorh,'Visible','on')
else
set(hpatch,'FaceVertexCData',vertexcolorh,'Visible','on')
end
end
vertexcolorh = []; % Move to next surface
figure(ImagingFigure)
hold off
rotate3d(PlotHandles.hFigure,'on')
if ~isempty(findobj(gcf,'Tag','MainAxes'))
if length(hpatch)>1 & length(imesh) == 1% suplots with only one surface to display (need to generalize this)
if ~strcmp(get(hpatch(4),'Visible'),'on')
set(gca,'View',get(findobj(gcf,'Tag','MainAxes'),'view'))
end
end
end
%dataplot_cb mesh_lighting_props
%------------------------------------ SubFunctions
function [] = DataOnSensors % Sensors are color-coded according to the data at each location
%dataplot_cb mapping_create
Visu = get(DATAPLOT,'Userdata');
previous = openfig('tessellation_window.fig','reuse');
% Handles to the spherical representation of the sensors
sph = findobj(previous,'Tag','SENSORS');
if isempty(sph) | isempty(get(previous,'userdata'))
if ~isempty(sph)
if strcmp(get(sph(1),'type'),'surface')
sph = dataplot_cb('see_sensors');
else
sph = dataplot_cb('see_sensors_markers');
end
else
TessWin = findobj(0,'Tag','tesselation_select');
hTessWin = guihandles(TessWin);
if get(hTessWin.Sensors3D,'Value')
sph = dataplot_cb('see_sensors');
else
sph = dataplot_cb('see_sensors_markers');
end
end
set(previous,'UserData',sph);
else
sph = get(previous,'Userdata'); % Need this and line above to keep sensor ordering correct - sph = findobj creates an arbitrary ordering otherwise
end
% Bring the time slider - call the mapping window if available
MAPPING = openfig('mapping.fig','reuse');
hMAPPING = guihandles(MAPPING);
% Good/Bad Channels
if ~isempty(find(Visu.Data.ChannelFlag == -1))
goodchannels = good_channel(Visu.Channel,Visu.Data.ChannelFlag,DataType{current}); % Discard bad channels
badchannels = setdiff(good_channel(Visu.Channel,ones(length(Visu.Channel),1),DataType{current}),goodchannels);
% goodchannels = setdiff(Visu.ChannelID{current},find(Visu.Data.ChannelFlag == -1));
%badchannels = intersect(Visu.ChannelID{current},find(Visu.Data.ChannelFlag == -1));
else
goodchannels = setdiff(Visu.ChannelID{current},find(Visu.Data.ChannelFlag == -1));
badchannels = [];
end
% Sampling rate
srate = abs(Visu.Data.Time(1)-Visu.Data.Time(2) ); % Assuming time is in sec.
ctime = round((get(hMAPPING.hDataplot.TimeSlider,'Value')/1000-Visu.Data.Time(1))/srate)+1;
axx = get(sph(1),'Parent');
data = Visu.Data.F(goodchannels,:); clear Visu
%set(axx,'Clim',[-max(abs(data(:))) max(abs(data(:)))]);
caxis(axx,[-max(abs(data(:))) max(abs(data(:)))]);
data = data(:,ctime);
C = [flipud(fliplr(hot(128)));hot(128)];
colormap(axx,C);
% Manipulation to get the color properly for each sensor sphere
set(sph,'visible','off')
TessWin = findobj(0,'Tag','tesselation_select');
hTessWin = guihandles(TessWin);
if get(hTessWin.Sensors3D,'Value')
zdata = get(sph(1),'Zdata');
for k=1:length(sph)
cdata = data(k)*ones(size(get(sph(k),'zdata')));
set(sph(k),'CData',cdata,'facecolor','flat','CDataMapping','scaled');
end
else
for k=1:length(sph)
set(sph(k),'Markerfacecolor',C( max([round(size(C,1)*(data(k)-min(caxis(axx)))/(max(caxis(axx))-min(caxis(axx)))),1]),:));
end
end
set(sph,'visible','on')
%--------------------------------------------------------------------------------------------------------------------
function [] = DataOnScalp % Interpolated Scalp Topography of Surface Data
%dataplot_cb mapping_create
TessSelect = findobj(0,'Tag','tesselation_select');
hSelect = guihandles(TessSelect);
ActiveTess = get(hSelect.removed,'String'); % Find the active scalp surface
% Check if the CORTEX keyword is present
for k = 1:length(ActiveTess)
flag(k) = ~isempty(findstr(lower(ActiveTess{k}),'head'));%'scalp'));
end
iScalp = find(flag);
iScalp = iScalp(1);
if isempty(iScalp)
h = msgbox('Please select a scalp envelope from the tessellated surfaces available.');
return
end
TessWin = openfig('tessellation_window.fig','reuse');
Visu = get(DATAPLOT,'Userdata');
if ~isempty(find(Visu.Data.ChannelFlag == -1))
goodchannels = good_channel(Visu.Channel,Visu.Data.ChannelFlag,DataType{current}); % Discard bad channels
badchannels = setdiff(good_channel(Visu.Channel,ones(length(Visu.Channel),1),DataType{current}),goodchannels);
% goodchannels = setdiff(Visu.ChannelID{current},find(Visu.Data.ChannelFlag == -1));% - min(Visu.ChannelID{current})+1; % Discard bad channels
%badchannels = intersect(Visu.ChannelID{current},find(Visu.Data.ChannelFlag == -1));% - min(Visu.ChannelID{current})+1;
else
goodchannels = Visu.ChannelID{current};
badchannels = [];
end
Scalp = findobj(TessWin,'Type','patch','Tag',ActiveTess{iScalp}); % Need to find a better way to identify scalp when multiple surfaces are present in tessellation_win
if isempty(Scalp), tesselation_select_done, return, end
if isempty(get(Scalp,'Userdata'))
figure(TessWin)
C = [flipud(fliplr(grayish(hot(128),.3)));grayish(hot(128),.3)];
colormap(C)
scalp.vertices = get(Scalp,'vertices');
% Channels Locations
j = 0;
for chan = goodchannels
j = j+1;
%Projection of sensors on surface
celec = ones(size(scalp.vertices,1),1)*Visu.Channel(chan).Loc(:,1)';
dist = norlig(scalp.vertices-celec);
[minn(j) imin(j)] = min(dist);
end
scalp.chanloc = [scalp.vertices(imin,1),scalp.vertices(imin,2),scalp.vertices(imin,3)];
% Sampling rate
scalp.srate = abs((Visu.Data.Time(2)-Visu.Data.Time(1)));
data = Visu.Data.F(goodchannels,:);
axx= get(Scalp,'Parent');
set(axx,'Clim',[-max(abs(data(:))) max(abs(data(:)))]);
set(Scalp,'Userdata',scalp)
else
scalp = get(Scalp,'Userdata');
end
% Get the data from the mapping window
MAPPING = openfig('mapping.fig','reuse');
hMAPPING = guihandles(MAPPING);
ctime = round((get(hMAPPING.hDataplot.TimeSlider,'Value')/1000-Visu.Data.Time(1))/scalp.srate)+1;
vertxcolor = interp_mail(scalp.vertices,scalp.chanloc,Visu.Data.F(goodchannels,ctime));
set(Scalp,'FaceVertexCData',vertxcolor,'facecolor','interp');
%--------------------------------------------------------------------------------------------------------------------
function [] = CorticalMap(TessWin,OPTIONS) % See cortical current maps interpolated on the proper cortical surface
% TessSelect = findobj(0,'Tag','tesselation_select');
% hSelect = guihandles(TessSelect);
% ActiveTess = get(hSelect.removed,'String'); % Find the active scalp surface
% iCortex = get(hSelect.removed,'Value'); % Find the active scalp surface
% if isempty(iCortex)
% h = msgbox('Please select a cortex surface from the tessellations available.');
% return
% end
%TessWin = findobj(0,'Tag','tessellation_window');
SlidesWin = findobj(0,'Tag','subplot_tess_window');
if ~isempty(SlidesWin)
TessWin = SlidesWin; % Subplot or movie window
end
% if isempty(TessWin)
% TessWin = openfig('tessellation_window.fig');
% dataplot_cb('loaddata',Users.CurrentData.StudyFile,bst_static_taskbar('GET','DATA')); % Refresh by reloading data
% end
% Check if loaded Results match this surface (ie number of sources = number of vertices of selected cortical surface)
Results = get(hSelect.ResultFiles,'Userdata');
nSources = size(Results.ImageGridAmp,1);
if nargin == 1 & ~OPTIONS.OrthoViews
if iscell(ActiveTess)
Cortex = findobj(get(TessWin,'CurrentAxes'),'Type','patch','Tag',ActiveTess{iCortex}); %CHEAT - need to be improved ?
else
Cortex = findobj(get(TessWin,'CurrentAxes'),'Type','patch','Tag',ActiveTess); %CHEAT - need to be improved ?
end
elseif nargin == 1 & get(hSelect.OrthoViews,'Value') == 1 % Orthogonal views are requested
Cortex = findobj(get(TessWin,'CurrentAxes'),'Type','patch','Tag',ActiveTess,'Visible','on'); %CHEAT - need to be improved ?
Cortex = Cortex(1);
else
Cortex = varargin{1}(1);
end
%set(Cortex,'Visible','off')
if ~isfield(get(Cortex,'Userdata'),'srate')
Visu = get(DATAPLOT,'Userdata');
figure(TessWin)
cortex = get(Cortex,'Userdata');
cortex.vertices = get(Cortex,'vertices');
% Sampling rate
cortex.srate = abs(Visu.Data.Time(2)-Visu.Data.Time(1));
data = Results.ImageGridAmp;
axx= get(Cortex,'Parent');
set(axx,'Clim',[-max(abs(data(:))) max(abs(data(:)))]);
set(Cortex,'Userdata',cortex)
% Result time window has priority over the original time window for the data - replace
Visu.Data.Time = Results.ImageGridTime; % Need to do better than that - CHEAT
try
Visu.Data.F = Visu.Data.F(:,Results.Time(1):Results.Time(end)); % Should work only the first time the result time series are called - otherwise do nothing
catch
Visu.Data.F = Visu.Data.F; % Do nothing
end
set(DATAPLOT,'Userdata',Visu);
M = max(abs(Results.ImageGridAmp(:)));
set(hSelect.Colorbar,'Userdata',M);
else
cortex = get(Cortex,'Userdata');
end
%dataplot_cb mapping_create
% Get the data from the mapping window
MAPPING = findobj(0,'Tag','mapping');
hMAPPING = guihandles(MAPPING);
if iscell(cortex)
ctime = round((get(hMAPPING.hDataplot.TimeSlider,'Value')/1000- Results.ImageGridTime(1))/cortex{1}.srate)+1;
else
ctime = round((get(hMAPPING.hDataplot.TimeSlider,'Value')/1000- Results.ImageGridTime(1))/cortex.srate)+1;
end
if ctime < 0, ctime =1; end
if get(hSelect.AbsoluteCurrent,'Value')
set(Cortex,'FaceVertexCData',abs(Results.ImageGridAmp(:,ctime)));
else
set(Cortex,'FaceVertexCData',Results.ImageGridAmp(:,ctime));
end
figure(TessWin)
hTessWin = guihandles(TessSelect);
if get(hTessWin.Normalize,'Value')
TessWin = findobj(0,'Tag','tessellation_window');
axx = findobj(TessWin,'tag','MainAxes');
set(axx,'ClimMode','auto')
M = max(abs(Results.ImageGridAmp(:,ctime)));
set(hSelect.Normalize,'Userdata',M)
SlidesWin = findobj(0,'Tag','subplot_tess_window');
if ~isempty(SlidesWin)
TessWin = SlidesWin; % Subplot or movie window
axx = findobj(TessWin,'type','axes');
set(axx,'ClimMode','auto')
end
else
if get(hTessWin.FreezeColormap,'Value')
data = get(Cortex,'FaceVertexCData');
FreezeTime = get(hTessWin.FreezeColormap,'Userdata');
M = max(abs(Results.ImageGridAmp(:,FreezeTime)));
else
M = max(abs(Results.ImageGridAmp(:)));
end
if get(hSelect.AbsoluteCurrent,'Value')
set(findobj(TessWin,'type','axes'),'ClimMode','manual','Clim',[0 M])
else
set(findobj(TessWin,'type','axes'),'ClimMode','manual','Clim',[-M M])
end
end
if isempty(get(hSelect.ColorMAP,'Userdata'))
%if ~isfield(cortex,'CDepth') & ~get(hSelect.MapCurvature,'Value')% No curvature mapping
if ~get(hSelect.MapCurvature,'Value')% No curvature mapping
if ~get(hSelect.AbsoluteCurrent,'Value')
C = [flipud(fliplr(grayish(hot(128),.3)));grayish(hot(128),.3)];
else
%C = [(grayish(hot(128),.3))];
load bst_cactivcmap
C = map;
end
else
[FVC,C]=catci(get(Cortex,'FaceVertexCData'),cortex.CDepth,M);
set(Cortex,'FaceVertexCData',FVC);
end
else % Truncated Colormap
if ~isfield(cortex,'CDepth') % No curvature mapping
C = get(hSelect.ColorMAP,'Userdata'); % Truncated Colormap
else
if get(hSelect.ZScoreThresholdApply,'Value') & get(hSelect.ZScore,'Value')
cThres = 100*get(hSelect.ZScoreThreshold,'Userdata');
else
cThres = str2num(get(hTessWin.TruncateFactor,'String'));
end
[FVC,C]=catci(get(Cortex,'FaceVertexCData'),cortex.CDepth,M,cThres);
if get(hSelect.OrthoViews,'Value') == 0
set(Cortex,'FaceVertexCData',FVC);
else
Cortex = findobj(TessWin,'Type','patch','Tag',ActiveTess,'Visible','on'); %CHEAT - need to be improved ?
set(Cortex,'FaceVertexCData',FVC);
end
end
end
colormap(C)
% % ------------------------
% case 'mapping_create' % Generates SLIDES and/or MOVIES
%
% Visu = get(DATAPLOT,'Userdata');
% MEG = get(findobj(DATAPLOT,'Tag','MEG'),'value');
% EEG = get(findobj(DATAPLOT,'Tag','EEG'),'value');
% OTHER = get(findobj(DATAPLOT,'Tag','OTHER'),'value');
% current = find([MEG,EEG,OTHER] == 1);
%
% mapping_win = findobj(0,'Tag','mapping');
% if isempty(mapping_win)
% open mapping.fig;
% mapping_win = findobj(0,'Tag','mapping');
% else
% return
% end
%
% RAW = findobj(mapping_win ,'Tag','raw');
% GRADIENT = findobj(mapping_win ,'Tag','gradient');
% MAGNETIC = findobj(mapping_win,'Tag','magnetic');
%
% if current == 1 % If MEG
% set([RAW,GRADIENT,MAGNETIC],'enable','on')
% set(RAW,'Value',1)
% else
% set([RAW,GRADIENT,MAGNETIC],'enable','off')
% end
%
% start = findobj(mapping_win,'Tag','start');
% step = findobj(mapping_win,'Tag','step');
% stop = findobj(mapping_win,'Tag','stop');
% current_time = findobj(mapping_win,'Tag','current_time');
% hDataplot.TimeSlider = findobj(mapping_win,'Tag','hDataplot.TimeSlider');
% go = findobj(mapping_win,'Tag','go');
% SPHERE = findobj(mapping_win,'Tag','sphere');
% fit = findobj(mapping_win,'Tag','fit');
% ROTATE = findobj(mapping_win,'Tag','rotate');
% fit = findobj(mapping_win,'Tag','fit');
% slides = findobj(mapping_win,'Tag','slides');
% single = findobj(mapping_win,'Tag','single');
% MOVIE = findobj(mapping_win,'Tag','movie');
%
% handles = guihandles(mapping_win);
%
% set(mapping_win,'Name',[get(mapping_win,'Name'), ' - ', modality{current}])
%
% set(start,'String', num2str(1000*Visu.Data.Time(1)))
% set(stop,'String', num2str(1000*Visu.Data.Time(end)))
% if length(Visu.Data.Time) > 1
% set(step,'String', num2str(1000*(Visu.Data.Time(2)-Visu.Data.Time(1)),2))
% else
% set(step,'String', 0)
% end
%
%
% data.TIME_MIN = str2num(get(findobj(DATAPLOT,'Tag','time_min'),'String'));
% data.TIME_MAX = str2num(get(findobj(DATAPLOT,'Tag','time_max'),'String'));
%
% if data.TIME_MAX > Visu.Data.Time(end) * 1000
% data.TIME_MAX = Visu.Data.Time(end) * 1000;
% set(findobj(DATAPLOT,'Tag','time_max'),'String',num2str(data.TIME_MAX,5))
% end
%
% if data.TIME_MIN < Visu.Data.Time(1) * 1000
% data.TIME_MIN = Visu.Data.Time(1) * 1000;
% set(findobj(DATAPLOT,'Tag','time_min'),'String',num2str(data.TIME_MIN,5))
% end
%
% if ~isempty(mapping_win)
% hDataplot.TimeSlider = findobj(mapping_win,'Tag','hDataplot.TimeSlider');
% if length(Visu.Data.Time) > 1
% set(hDataplot.TimeSlider,'enable','on')
% set(hDataplot.TimeSlider,'Min',data.TIME_MIN,'Max',data.TIME_MAX)
% else
% set(hDataplot.TimeSlider,'Min',data.TIME_MIN,'Max',2*data.TIME_MAX)
% set(hDataplot.TimeSlider,'enable','off')
% end
%
% end
% set(hDataplot.TimeSlider,'Value',data.TIME_MIN,'enable','on')
%
% set(current_time,'String',num2str(get(hDataplot.TimeSlider,'Value'),4))
% set(single,'Value',1)
% mutincomp([single,slides,MOVIE])
% set(SPHERE,'Value',1)
% mutincomp([SPHERE,fit,handles.scalp])
%
% % Draw time cursor
% figsingle = findobj(0,'Tag',['waves_single',int2str(current)]); % Figure of overlaping plots for the current modality
% if ~isempty(figsingle)
% figure(figsingle)
% ctime = get(hDataplot.TimeSlider,'Value');
% cursor = line([ctime ctime],get(gca,'Ylim'));
% set(cursor,'color','r','linewidth',3,'Tag','cursor','erasemode','Xor')
% end
%
% fignplot = findobj(0,'Tag','nplot_win'); % superimposed n-plots
% if ~isempty(fignplot)
% figure(fignplot)
% ctime = get(hDataplot.TimeSlider,'Value');
% cursor = line([ctime ctime],get(gca,'Ylim'));
% set(cursor,'color','r','linewidth',3,'Tag','cursor','erasemode','Xor')
% end
% %-----------------------------------------------------------------------------------------------
function [mixedRGB] = mix_curvature_current(varargin)
% Various colormaps ---------------------------------------
data.curDensityMapArray{1} = [zeros(17,2) (0.5:0.03:1)'
zeros(29,1) (0.03:0.034:1)' ones(29,1)
(0.3:0.3:1)' ones(3,2)
1 1 1
ones(26,1) (1:-0.039:0)' (1:-0.039:0)'
(1:-0.021:0.5)' zeros(24,2)];
data.curDensityMapArray{2} = [zeros(17,2) (0.5:0.03:1)' %17
zeros(29,1) (0.03:0.034:1)' ones(29,1) %29
(0.3:0.3:1)' ones(3,2) %3
1 1 1
ones(50,2) (1:-0.0202:0)']; %50
tmp = hot(140);%flipdim(hot(100),1);
data.curDensityMapArray{3} = tmp(1:100,:);%flipdim(hot(100),1);
% Various colormaps ------------------DONE----------------
[hf,hs,hl] = view_surface('test',Faces,tmpVertices,Results.ImageGridAmp(:,1));
delete(findobj(gcf,'type','light'))
%view(-10,0),
material dull, camlight
set(hs,'backfacelighting','unlit')
hold on
%set(hf,'color','k')
colormap flag; colormenu
flag = get(gcf,'colormap');
flag = flag(:,[3 2 1]);
colormap(flag)
timee = 1:timestep:size(Results.ImageGridAmp,2)-corrsize/2;%
corr = zeros(size(Results.ImageGridAmp(patchh,:),1),length(timee));
tt = 0;
tmp = Results.ImageGridAmp(patchh,timee);
corr = abs(tmp)/max(tmp(:));%abs(activ(:,1+corrsize/2));
corr(corr<cd_threshold) = 0; %Show current density that is at least 'display_threshold'% the maximum absolute value
limitValue = max(abs(corr(:)));
%Apply alpha blending
%create anatomy and current density colormaps
if(anatomy==1) anatomy=0.99; end
scale = linspace(.4,.6,100);%anatomy:(1-anatomy)/99:1;
data.anatMap = repmat(scale',1,3);
if(viewmode == 1) %if real values
data.curDensityMap = data.curDensityMapArray{colormap_selected};
else %if abs values
data.curDensityMap = data.curDensityMapArray{colormap_selected+2};
end
colormap(data.curDensityMap);
%Get RGB info for anatomy
min_curv = min(curvature);
max_curv = max(curvature);
data.index_anatomy = (floor((curvature-min_curv)*99/(max_curv-min_curv)+1));
data.anatRGB = data.anatMap(data.index_anatomy,:);
for t = 1:length(timee)
% tmp = NaN * Results.ImageGridAmp(:,t);
% tmp(patchh) = corr(:,t);
%
if(viewmode == 1)%if real values
index_density = floor((corr(:,t)+limitValue)*99/(2*limitValue)+1);
data.curDensityMap = data.curDensityMapArray{colormap_selected};
else %if abs values
corr = abs(corr);
index_density = floor((corr(:,t)-limitValue)*99/limitValue+100);
data.curDensityMap = data.curDensityMapArray{colormap_selected+2};
end
data.curDensityRGB = data.curDensityMap(index_density,:);
data.toBlend = find(corr(:,t) ~= 0);
mixedRGB = data.anatRGB;
toBlend = data.toBlend;
mixedRGB(patchh(toBlend),:) = cd_transparency*data.anatRGB(toBlend,:)+(1-cd_transparency)*data.curDensityRGB(toBlend,:);
set(hs,'FaceVertexCData',mixedRGB);
drawnow
end
Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Oct 12 12:05:14 2004
Cross-Directory links are: ON