Description of vb_load

0001 function [B, Model, nGact, nGall, nGGall, Cov, Wact, Wbck, ...
0002              bsnorm, ix_act, ix_bck, vb_parm] ...
0003          = vb_load_result(proj_root, resultfile, ix_area)
0004 % Load VB estimation result and necessary data for current reconstruction
0005 %   [B, Model, nGact, nGall, nGGall, Cov, Wact, Wbck, ...
0006 %       bsnorm, ix_act, ix_bck, vb_parm] ...
0007 %   = vb_load_result(proj_root, resultfile, ix_area)
0008 % --- Input
0009 % proj_root
0010 % resultfile : result file obtained by VB estimation
0011 % ix_area : Position index to calculate estimated current
0012 %    If 'ix_area' is empty or not given,
0013 %    currents in the focal region are calculated
0014 %
0015 % 2006-08-31 M. Sato
0016 % 2006-11-11 M. Sato
0017 % If 'ix_area' is empty, original focal area is used
0018 % instead of extended focal area by spatial smoothing
0019 %
0020 % Copyright (C) 2011, ATR All Rights Reserved.
0021 % License : New BSD License(see VBMEG_LICENSE.txt)
0022 
0023 
0024 if ~isempty(proj_root)
0025   resultfile = [proj_root filesep resultfile];
0026 end  
0027 
0028 load(resultfile, 'bayes_parm','Model','Cov','vb_parm');
0029 
0030 if ~isempty(proj_root);
0031   bayes_parm_abs = vb_parm_absolute_path(proj_root, bayes_parm);
0032 else
0033   bayes_parm_abs = bayes_parm;
0034 end
0035 
0036 if isfield(bayes_parm_abs,'remove_crossed_area'),
0037     remove_crossed_area = bayes_parm_abs.remove_crossed_area;
0038 else
0039     remove_crossed_area = OFF;
0040 end
0041 
0042 % MEG data preparation
0043 [B] = vb_megdata_preparation(bayes_parm_abs);
0044 
0045 % Preparation of lead fields
0046 [nGall, bsnorm, ix_global, ix_global_ex, Wall, Lbck] = ...
0047     vb_global_leadfield_preparation(bayes_parm_abs);
0048 [nGact, ix_focal, ix_focal_ex, Wact, Lact]  =  ...
0049     vb_focal_leadfield_preparation(bayes_parm_abs, bsnorm);
0050 
0051 if remove_crossed_area == ON,
0052     % Remove focal area from global area
0053     switch bayes_parm_abs.forward_model
0054     case 'focal+global',
0055         [nGall, bstmp, ix_global, ix_global_ex, Wall, Lbck] = ...
0056             vb_global_leadfield_preparation(bayes_parm_abs, bsnorm, ON );
0057     case 'focal',
0058         ix_global_ex = [];
0059     end
0060 end
0061 
0062 %%%% Find active current area
0063 
0064 % Current variance
0065 A_inv = sum(Model.a , 2);    % sum over all time window
0066 Nvact = vb_parm.Nvact;        % # of active vertex
0067 
0068 % Multiply smoothing filter
0069 if bayes_parm_abs.variance_orientation == ON
0070     A_inv = reshape(A_inv, [Nvact Lact]);
0071     W = sum( Wact * A_inv ,2);
0072 else
0073     W = Wact*A_inv;
0074 end
0075 
0076 % Find nonzero component of (Wact * A_inv) : active current region
0077 ix_act_nz   = find( W > 0 );
0078 ix_focal_nz = ix_focal_ex(ix_act_nz);
0079 
0080 % Area in which current is calculated
0081 if isempty(ix_area),
0082     ix_area = ix_focal;        % original focal area
0083 %    ix_area = ix_focal_nz;    % extended focal area by spatial smoothing
0084 end
0085 
0086 % Select vertex index 'ix_area' within the active current region
0087 [j_tmp, ix_area_nz] = vb_index2indexsub(ix_area, ix_focal_nz);
0088 
0089 %%%% Area information where the currents are reconstructed.
0090 
0091 % Convert vertex index to index within the active current
0092 % Indices outside 'ix_focal_ex' are removed from current indices 'ix_area'
0093 
0094 % ix_act : current indices w.r.t all vertices
0095 % j_act  : current indices 'ix_act' w.r.t ix_focal_ex
0096 [j_act, ix_act] = vb_index2indexsub(ix_area_nz, ix_focal_ex);
0097 
0098 % Global(background) area is independent of active current region
0099 if remove_crossed_area == OFF
0100     % select active region
0101     [k_bck, ix_bck] = vb_index2indexsub(ix_area_nz, ix_global_ex);
0102 else
0103     [k_bck, ix_bck] = vb_index2indexsub(ix_area, ix_global_ex);
0104 end
0105 
0106 % # of focal window
0107 Lact     = vb_parm.Norient;
0108 Lact_var = vb_parm.Norient_var;
0109 Njact     = vb_parm.Njact;            % # of active current parameters
0110 Nvact_area = length(ix_act);        % # of active vertex in area
0111 Njact_area = Nvact_area * Lact;     % # of current parameters in area
0112 
0113 % # of global window
0114 Lbck  = vb_parm.Norient_all;
0115 Nvall = length(ix_global);         % # of dipoles in the background
0116 Nvall_area = length(ix_bck);       % # of background points in area
0117 Njall_area = Nvall_area* Lbck ;    % # of current parameters in area
0118 
0119 vb_parm.Njact_area = Njact_area;
0120 vb_parm.Njall_area = Njall_area;
0121 
0122 fprintf('# of active current vertex = %d\n', Nvact_area)
0123 fprintf('# of background vertex = %d\n', Nvall_area)
0124 
0125 Nsession = vb_parm.Nsession;    % Number of sessions
0126 
0127 %%  Gb * Gb'
0128 for n = 1 : Nsession
0129      nGGall{n} = nGall{n}*nGall{n}';    
0130 end
0131 
0132 %% Restrict smoothing filter 'W' to estimation points
0133 Wact = Wact(j_act,:);    % Nvact_area x Nvact
0134 Wbck = Wall(k_bck,:);   % Nvall_area x Nvall
0135 
0136 %% Gb*Wb' restricted to ix_area
0137 %% --- For faster computation ---
0138 for n=1:Nsession
0139     Gb   = nGall{n};   % Lead field of background
0140     Nch  = size(Gb,1); % # of sensor
0141     
0142     GbW  = reshape(Gb,  [Nch*Lbck Nvall]) * Wbck'; % [Nch*Lbck x Nvall_area]
0143     GbW  = reshape(GbW, [Nch Lbck*Nvall_area]);    % [Nch x Lbck*Nvall_area]
0144     
0145     nGall{n} = GbW;
0146 end
0147 
0148 %    switch Lbck
0149 %    case 1,
0150 %        GbW  =  Gb*Wbck';
0151 %    case 2,
0152 %        GbW  = [Gb(:,1:Nvall)*Wbck', Gb(:,Nvall+1:2*Nvall)*Wbck'];
0153 %    case 3,
0154 %        GbW  = [Gb(:,1:Nvall)*Wbck', Gb(:,Nvall+1:2*Nvall)*Wbck', ...
0155 %                Gb(:,2*Nvall+1:3*Nvall)*Wbck'];
0156 %    end
0157 
0158 clear Gb;
vb_load_result

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
