Localize scale

best/cmem/clusters/be_spatial_priorw.m

@@ -31,8 +31,7 @@
 %    You should have received a copy of the GNU General Public License
 %    along with BEst. If not, see <http://www.gnu.org/licenses/>.
 % -------------------------------------------------------------------------
-
-
+
 % Default options settings
 Def_OPTIONS.solver = struct('spatial_smoothing', 0.6);
 
@@ -50,13 +49,13 @@
 rho = OPTIONS.solver.spatial_smoothing; % scalar that weight the adjacency matrix 
 W   = speye(nb_vertices);
 
+% If we don't smooth, or if neighbor is emtpy, nothing to do.
+if (rho == 0 || isempty(neighbors))
+    return;
+end
 
 % Add comment to result
-if rho && ~isempty(neighbors)
-    OPTIONS.automatic.Comment = [OPTIONS.automatic.Comment ' | smooth=' num2str(rho)];
-else
-    return
-end
+OPTIONS.automatic.Comment = [OPTIONS.automatic.Comment ' | smooth=' num2str(rho)];
 
 % Preallocation of the sparse matrix
 A  = neighbors - spdiags(sum(neighbors,2),0,nb_vertices,nb_vertices);

best/cmem/preprocess/be_normalize_and_units.m

@@ -56,7 +56,8 @@
     %% ===== Compute Minimum Norm Solution ==== %% 
     % we compute MNE (using l-curve for nirs or depth-weighted version)
 
-    [obj, OPTIONS] = be_main_mne(obj, OPTIONS);
+    OBJ_FUS = be_fusion_of_modalities(obj, OPTIONS, 0);
+    [~, OPTIONS] = be_main_mne(OBJ_FUS, OPTIONS);
 
 end
 

best/cmem/solver/be_cmem_solver.m

@@ -89,15 +89,14 @@
 % Convert to average reference (only for EEG / iEEG)
 [OPTIONS]       = be_avg_reference(OPTIONS);
 
-%% ===== Pre-whitening of the data ==== %%
-% it uses empty-room data if available
-% if PlOS one : nothing is done here
-% [OPTIONS] = be_prewhite(OPTIONS);
-
 %% ===== Pre-process the leadfield(s) ==== %% 
 % we keep leadfields of interest; we compute svd of normalized leadfields
 [OPTIONS, obj] = be_main_leadfields(obj, OPTIONS);
 
+%% ===== pre-processing for spatial smoothing (Green mat. square) ===== %%
+% matrix W'W from the Henson paper
+[OPTIONS, obj.GreenM2] = be_spatial_priorw( OPTIONS, obj.VertConn);
+
 %% ===== Apply temporal data window  ===== %%
 % check for a time segment to be localized
 [OPTIONS] = be_apply_window(OPTIONS, []);
@@ -107,19 +106,19 @@
 % we normalize the data and the leadfields
 [OPTIONS, obj] = be_normalize_and_units(obj, OPTIONS);
 
-%% ===== Clusterize cortical surface ===== %%
-[OPTIONS, obj]  = be_main_clustering(obj, OPTIONS);
-
-%% ===== pre-processing for spatial smoothing (Green mat. square) ===== %%
-% matrix W'W from the Henson paper
-[OPTIONS, obj.GreenM2] = be_spatial_priorw( OPTIONS, obj.VertConn);
-
 %% ===== Noise estimation ===== %%   
 [OPTIONS, obj] = be_main_data_preprocessing(obj, OPTIONS);
 
+%% ===== Clusterize cortical surface ===== %%
+[OPTIONS, obj]  = be_main_clustering(obj, OPTIONS);
+
 %% ===== Fuse modalities ===== %%   
 obj = be_fusion_of_modalities(obj, OPTIONS);
-
+
+%% ===== Set Alpha ===== %%
+% Set Alpha value for each cluster
+[OPTIONS, obj]  = be_main_alpha(obj, OPTIONS);
+
 %% ===== Solve the MEM ===== %%
 [obj.ImageGridAmp, OPTIONS] = be_launch_mem(obj, OPTIONS);
 if OPTIONS.optional.display
@@ -130,7 +129,7 @@
 [obj, OPTIONS] = be_unormalize_and_units(obj, OPTIONS);
 
 %% ===== Inverse temporal data window  ===== %%
-[OPTIONS, obj] = be_apply_window( OPTIONS, obj );
+[OPTIONS, obj] = be_apply_window(OPTIONS, obj);
 
 %% ===== Prepare output  ===== %%
 Results = be_template('resultsmat');

best/main/be_main_alpha.m

@@ -0,0 +1,73 @@
+function [OPTIONS, obj] = be_main_alpha(obj, OPTIONS)
+% BE_MAIN_ALPHA initialize the alpha value for each cluster
+%
+% Inputs:
+% -------
+%
+%	obj			:	MEM obj structure
+%   OPTIONS     :   structure (see bst_sourceimaging.m)
+%
+%
+% Outputs:
+% --------
+%
+%   OPTIONS     :   Updated options fields
+%	obj			:	Updated structure
+%
+%% ==============================================   
+% Copyright (C) 2011 - LATIS Team
+%
+%  Authors: LATIS team, 2011
+%
+%% ==============================================
+% License 
+%
+% BEst is free software: you can redistribute it and/or modify
+%    it under the terms of the GNU General Public License as published by
+%    the Free Software Foundation, either version 3 of the License, or
+%    (at your option) any later version.
+%
+%    BEst is distributed in the hope that it will be useful,
+%    but WITHOUT ANY WARRANTY; without even the implied warranty of
+%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%    GNU General Public License for more details.
+%
+%    You should have received a copy of the GNU General Public License
+%    along with BEst. If not, see <http://www.gnu.org/licenses/>.
+% -------------------------------------------------------------------------   
+
+if ~isfield(OPTIONS.optional.clustering, 'initial_alpha')
+
+    %% ===== Double to single precision  ===== %%
+    [OPTIONS] = be_switch_precision(OPTIONS, 'single');
+
+    if OPTIONS.model.alpha_method < 6
+        [ALPHA, CLS, OPTIONS] = be_scores2alpha(obj.SCR, obj.CLS, OPTIONS);
+    else % We compute the score using MNE
+        [ALPHA, CLS, OPTIONS] = be_mne2alpha(obj , obj.CLS, OPTIONS);
+    end
+
+    %% ===== Single to double precision  ===== %%
+    [OPTIONS] = be_switch_precision(OPTIONS, 'double');
+
+
+else
+    CLS = obj.CLS;
+    if strcmp( OPTIONS.mandatory.pipeline, 'wMEM' )
+        ALPHA = OPTIONS.optional.clustering.initial_alpha * ones(1,size(OPTIONS.automatic.Modality(1).selected_jk, 2));
+    else
+        ALPHA = OPTIONS.optional.clustering.initial_alpha * ones(1,size(OPTIONS.automatic.Modality(1).data, 2));
+    end
+end
+
+% the final clusters (CLS) and alpha's (ALPHA)
+obj.CLS   = CLS;
+obj.ALPHA = ALPHA;
+
+end
+
+
+
+
+
+

best/main/be_main_clustering.m

@@ -39,50 +39,35 @@
 
 if ~isfield(OPTIONS.optional.clustering, 'initial_alpha')
     % Sources prescoring - MSP (ref. Mattout et al. 2006) and clustering
-    % following different strategies depending of the case cmem, wmem or
-    % rmem
-
+
     %% ===== Double to single precision  ===== %%
     [OPTIONS] = be_switch_precision(OPTIONS, 'single');
 
     switch OPTIONS.mandatory.pipeline
         case 'cMEM'
             [CLS, SCR, OPTIONS] = be_cmem_clusterize_multim(obj, OPTIONS); 
-            [ALPHA, CLS, OPTIONS]   = be_scores2alpha(SCR, CLS, OPTIONS);
         case 'wMEM'
             [CLS, SCR, OPTIONS] = be_wmem_clusterize_multim(obj, OPTIONS);
-
-            if OPTIONS.model.alpha_method < 6
-                [ALPHA, CLS, OPTIONS] = be_scores2alpha(SCR, CLS, OPTIONS);
-            else % We compute the score using MNE
-                [ALPHA, CLS, OPTIONS] = be_mne2alpha(obj , CLS, OPTIONS);
-            end
-
         case 'rMEM'
             [CLS, SCR, OPTIONS] = be_rmem_clusterize_multim(obj, OPTIONS);
-            [ALPHA, CLS, OPTIONS]   = be_scores2alpha(SCR, CLS, OPTIONS);
-
     end   
 
     %% ===== Single to double precision  ===== %%
     [OPTIONS] = be_switch_precision(OPTIONS, 'double');
 
 
 elseif strcmp( OPTIONS.mandatory.pipeline, 'wMEM' )
-    ALPHA = OPTIONS.optional.clustering.initial_alpha * ones(1,size(OPTIONS.automatic.Modality(1).selected_jk, 2));
     CLS   = OPTIONS.optional.clustering.clusters * ones(1,size(OPTIONS.automatic.Modality(1).selected_jk, 2));
     SCR   = [];
 
 else
-    ALPHA = OPTIONS.optional.clustering.initial_alpha * ones(1,size(OPTIONS.automatic.Modality(1).data, 2));
     CLS   = OPTIONS.optional.clustering.clusters * ones(1,size(OPTIONS.automatic.Modality(1).data, 2));
     SCR   = [];
 end
 
 % the final scores (SCR), clusters (CLS) and alpha's (ALPHA)
 obj.SCR   = SCR;
 obj.CLS   = CLS;
-obj.ALPHA = ALPHA;
 
 end
 

best/misc/be_fusion_of_modalities.m

@@ -55,22 +55,40 @@
 data = [];
 data_normalized = [];
 
+scaling_data = [];
+scaling_data_normalized = [];
+
 for ii=1:length(OPTIONS.mandatory.DataTypes)
+
     if isfield(obj, 'data') % wavelet
         data_mod = obj.data{ii};
     else % Time-series
         data_mod = OPTIONS.automatic.Modality(ii).data;
     end
-    
+
     data = vertcat(data, data_mod);
     data_normalized = vertcat(data_normalized, bsxfun(@rdivide, data_mod, sqrt(sum(data_mod.^2, 1))));
+
+    % Concatenate scaling data if present
+    if isfield(obj,'scaling_data')
+        data_mod = obj.scaling_data{ii};
+
+        scaling_data   = vertcat(scaling_data,  data_mod);
+        scaling_data_normalized = vertcat(scaling_data_normalized, bsxfun(@rdivide, data_mod, sqrt(sum(data_mod.^2, 1))));
+    end
 end
 % remove nan from normalized data
 data_normalized(isnan(data_normalized)) = 0;
+scaling_data_normalized(isnan(scaling_data_normalized)) = 0;
 
 obj.data = data;
 obj.data_normalized = data_normalized;
 
+if isfield(obj,'scaling_data')
+    obj.scaling_data = scaling_data;
+    obj.scaling_data_normalized = scaling_data_normalized;
+end
+
 % Concatenate idata(complex data) if present
 if isfield(OPTIONS.automatic.Modality(1),'idata')
     obj.idata   = vertcat(OPTIONS.automatic.Modality.idata);

best/misc/be_fusion_of_samples.m

@@ -69,11 +69,10 @@
 
 else % only one modality, we keep the ordering related to the power 
     nboxes = size(OPTIONS.automatic.Modality(1).selected_jk,2);
-    OPTIONS.automatic.selected_samples = ...
-        [OPTIONS.automatic.Modality(1).selected_jk ; ones(1,nboxes)];
+    OPTIONS.automatic.selected_samples = [OPTIONS.automatic.Modality(1).selected_jk ; ones(1,nboxes)];
     % we may clear the OPTIONS.automatic.Modality(1).selected_jk
     % we code the modality in the last line of the table (1 here)
 end
 
 
-return
+end

best/mne/main/be_main_mne.m

@@ -12,21 +12,18 @@
 %       -   obj
 
     if nargin < 3 || isempty(method)
-        if OPTIONS.model.depth_weigth_MNE > 0 || any(strcmp( OPTIONS.mandatory.DataTypes,'NIRS')) 
+        if OPTIONS.model.depth_weigth_MNE > 0 || any(strcmp(OPTIONS.mandatory.DataTypes, 'NIRS')) 
             method = 'mne_lcurve';
         else
             method = 'mne';
         end
     end
 
-    % apply the fusion of modalities
-    OBJ_FUS = be_fusion_of_modalities(obj, OPTIONS, 0);
-
     switch(method)
         case 'mne_lcurve'
-            Kernel   =   be_jmne_lcurve(OBJ_FUS, OPTIONS, struct('hfig',obj.hfig, 'hfigtab',obj.hfigtab)); 
+            Kernel   =   be_jmne_lcurve(obj, OPTIONS, struct('hfig',obj.hfig, 'hfigtab',obj.hfigtab)); 
         case 'mne'
-            Kernel   =   be_jmne(OBJ_FUS, OPTIONS);
+            Kernel   =   be_jmne(obj, OPTIONS);
     end
 
     for ii  =   1 : numel(OPTIONS.mandatory.DataTypes)

best/mne/solver/be_cmne_solver.m

@@ -70,12 +70,13 @@
 % and the leadfields
 [OPTIONS, obj] = be_normalize_and_units(obj, OPTIONS);
 
+%% ===== Fuse modalities ===== %%   
+obj = be_fusion_of_modalities(obj, OPTIONS, 0);
+
 %% ===== Compute Minimum Norm Solution ==== %% 
 [obj, OPTIONS] = be_main_mne(obj, OPTIONS, OPTIONS.solver.mne_method);
 
 %% ===== Results  ===== %%
-OBJ_FUS = be_fusion_of_modalities(obj, OPTIONS, 0);
-
 obj.ImageGridAmp = OPTIONS.automatic.Modality.MneKernel * OBJ_FUS.data;
 [OPTIONS, obj]   = be_apply_window( OPTIONS, obj );
 

best/mne/solver/be_jmne_lcurve.m

@@ -155,8 +155,8 @@
         plot(ax, Prior, Fit,'b.');
         plot(ax, Prior(Index), Fit(Index),'ro');
         hold(ax,'off');
-        xlabel('Norm |WJ|');
-        ylabel('Residual |M-GJ|');
+        xlabel(ax, 'Norm |WJ|');
+        ylabel(ax, 'Residual |M-GJ|');
     end
 
 end

best/rmem/solver/be_rmem_solver.m

@@ -127,9 +127,15 @@
     [OPTIONS2, obj] = be_main_data_preprocessing(obj, OPTIONS2);
 
     %% ===== Clusterize cortical surface ===== %%
-
     [OPTIONS2, obj] = be_main_clustering(obj, OPTIONS2);
+
+    %% ===== Fuse modalities ===== %%   
+    obj = be_fusion_of_modalities(obj, OPTIONS);
 
+    %% ===== Set Alpha ===== %%
+    % Set Alpha value for each cluster
+    [OPTIONS, obj] = be_main_alpha(obj, OPTIONS);
+
     %% ===== pre-processing for spatial smoothing (Green mat. square) ===== %%
     % matrix W'W from the Henson paper
     [OPTIONS2, obj.GreenM2] = be_spatial_priorw( OPTIONS2, obj.VertConn);
@@ -166,7 +172,7 @@
     pause(.5); close(obj.hmem);
 end
 
-if ~OPTIONS.automatic.stand_alone | OPTIONS.automatic.process                                 
+if ~OPTIONS.automatic.stand_alone || OPTIONS.automatic.process                                 
 
  	% If added to a 'default_study' node: need to update results links12 30
 	OPTIONS.automatic           =   OPTIONS2.automatic;