vb_load_result_z

 Load VB estimation result and necessary data for current reconstruction
 [B, Model, nGact, nGall, nGGall, Cov, bsnorm, ix_act, ix_bck, ...
     vb_parm, Wact, Wbck, ix_act_ex, ix_bck_ex] ...
         = vb_load_result(proj_root, resultfile, ix_area)
 --- Input
 proj_root
 resultfile : result file obtained by VB estimation
 ix_area : Position index to calculate estimated current
    If 'ix_area' is empty or not given, 
    currents in the focal region are calculated

 2006-08-31 M. Sato
 2006-11-11 M. Sato
 If 'ix_area' is empty, original focal area is used
 instead of extended focal area by spatial smoothing

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function [B, Model, nGact, nGall, nGGall, Cov, bsnorm, ix_act, ix_bck, ...
0002              vb_parm, Wact, Wbck, ix_act_ex, ix_bck_ex] ...
0003          = vb_load_result_z(proj_root, resultfile, ix_area)
0004 % Load VB estimation result and necessary data for current reconstruction
0005 % [B, Model, nGact, nGall, nGGall, Cov, bsnorm, ix_act, ix_bck, ...
0006 %     vb_parm, Wact, Wbck, ix_act_ex, ix_bck_ex] ...
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 
0046 % Preparation of lead fields
0047 [nGall, bsnorm, ix_global, ix_global_ex, Wbck, Lbck] = ...
0048     vb_global_leadfield_preparation(bayes_parm_abs);
0049 [nGact, ix_focal, ix_focal_ex, Wact, Lact]  =  ...
0050     vb_focal_leadfield_preparation(bayes_parm_abs, bsnorm);
0051 
0052 if remove_crossed_area == ON,
0053     % Remove focal area from global area
0054     switch bayes_parm_abs.forward_model
0055     case 'focal+global',
0056         [nGall, bstmp, ix_global, ix_global_ex, Wbck, Lbck] = ...
0057             vb_global_leadfield_preparation(bayes_parm_abs, bsnorm, ON );
0058     case 'focal',
0059         ix_global = [];
0060         ix_global_ex = [];
0061     end
0062 end
0063 
0064 %%%% Find active current area
0065 
0066 % Current variance
0067 A_inv = sum(Model.a , 2);    % sum over all time window
0068 Nvact = vb_parm.Nvact;        % # of active vertex
0069 Nact  = length(A_inv);
0070 L     = Nact/Nvact;
0071 Lact  = vb_parm.Norient;
0072 
0073 if bayes_parm_abs.variance_orientation == ON
0074     A_inv = reshape(A_inv, [Nvact  L]);
0075     A_inv = sum( A_inv ,2);
0076 end
0077 
0078 % Find nonzero component of (A_inv) : active current region
0079 ix_act_nz   = find( A_inv > 0 );
0080 ix_focal_nz = ix_focal(ix_act_nz);
0081 
0082 % Area in which current is calculated
0083 if isempty(ix_area),
0084     ix_area = ix_focal;        % original focal area
0085 %    ix_area = ix_focal_nz;    % extended focal area by spatial smoothing
0086 end
0087 
0088 % Select vertex index 'ix_area' within the active current region
0089 [j_tmp, ix_area_nz] = vb_index2indexsub(ix_area, ix_focal_nz);
0090 
0091 %%%% Area information where the currents are reconstructed.
0092 
0093 % Convert vertex index to index within the active current
0094 % Indices outside 'ix_focal' are removed from current indices 'ix_area'
0095 
0096 % ix_act : current indices w.r.t all vertices
0097 % j_act  : current indices 'ix_act' w.r.t ix_focal
0098 [j_act, ix_act] = vb_index2indexsub(ix_area_nz, ix_focal);
0099 
0100 % Global(background) area is independent of active current region
0101 if remove_crossed_area == OFF
0102     % select active region
0103     [k_bck, ix_bck] = vb_index2indexsub(ix_area_nz, ix_global);
0104 else
0105     [k_bck, ix_bck] = vb_index2indexsub(ix_area, ix_global);
0106 end
0107 
0108 %Njact     = vb_parm.Njact;            % # of active current parameters
0109 
0110 % # of focal window
0111 Nvact_area = length(ix_act);        % # of active vertex in area
0112 Njact_area = Nvact_area * Lact;     % # of current parameters in area
0113 
0114 % # of global window
0115 Lbck  = vb_parm.Norient_all;
0116 Nvall = length(ix_global);         % # of dipoles in the background
0117 Nvall_area = length(ix_bck);       % # of background points in area
0118 Njall_area = Nvall_area* Lbck ;    % # of current parameters in area
0119 
0120 % Extract active current from a_inv
0121 a_inv = Model.a;
0122 Nwin  = size(a_inv,2);
0123 
0124 if bayes_parm_abs.variance_orientation == ON
0125     jj_act = repmat(j_act(:), [1 L]) ...
0126            + repmat((0:(L-1))*Nvact, [Nvact_area 1]);
0127 else
0128     jj_act = j_act;
0129 end
0130 
0131 Model.a = a_inv(jj_act(:),:);
0132 
0133 vb_parm.Njact_area = Njact_area;
0134 vb_parm.Njall_area = Njall_area;
0135 
0136 Nsession = vb_parm.Nsession;    % Number of sessions
0137 
0138 % Extract active current from nGact
0139 jj_act = repmat(j_act(:), [1 Lact]) ...
0140        + repmat((0:(Lact-1))*Nvact, [Nvact_area 1]);
0141 kk_bck = repmat(k_bck(:), [1 Lbck]) ...
0142        + repmat((0:(Lbck-1))*Nvall, [Nvall_area 1]);
0143 
0144 for n = 1 : Nsession
0145     G  = nGact{n};
0146     Gb = nGall{n};
0147 
0148     nGact{n}  =  G(:, jj_act(:));
0149     nGall{n}  = Gb(:, kk_bck(:));
0150      nGGall{n} = Gb * Gb';    
0151 end
0152 
0153 %% Restrict smoothing filter 'W' to estimation points
0154 Wact = Wact(:,j_act);    
0155 Wbck = Wbck(:,k_bck);   
0156 
0157 j_act_ex = find( sum(Wact,2) > 0);
0158 k_bck_ex = find( sum(Wbck,2) > 0);
0159 
0160 % ix_act_ex : current indices w.r.t all vertices
0161 % j_act_ex  : current indices w.r.t ix_focal_ex
0162 
0163 ix_act_ex = ix_focal_ex(j_act_ex);
0164 ix_bck_ex = ix_global_ex(k_bck_ex);
0165 
0166 Wact = Wact(j_act_ex,:);
0167 Wbck = Wbck(k_bck_ex,:);
0168 
0169 fprintf('# of active vertex   = %d\n', Nvact_area)
0170 fprintf('# of active current  = %d\n', Njact_area)
0171 
0172 fprintf('# of background vertex  = %d\n', Nvall_area)
0173 fprintf('# of background current = %d\n', Njall_area)
0174 
0175 fprintf('# of extended active vertex   = %d\n', length(ix_act_ex))
0176 fprintf('# of extended background vertex   = %d\n', length(ix_bck_ex))
0177 
0178 % ix_act : Vertex index corresponding to active current Zact
0179 % ix_bck : Vertex index corresponding to background current Zbck
0180 % ix_act_ex : Vertex index corresponding to active current Jact
0181 % ix_bck_ex : Vertex index corresponding to background current Jbck
0182 % Wact   : Spatial smoothing matrix of focal window
0183 % Wbck   : Spatial smoothing matrix of global window
0184 % Jact   = Wact * Zact
0185 % Jbck   = Wbck * Zbck
0186 
0187 % [Njact, Tall, Ntry] = size(Zact);
0188 % Njact_ex = size(Wact,1);
0189 %
0190 % Jact  = Wact * reshape(Zact,  [Njact/Lact, Lact*Tall*Ntry]);
0191 % Jact  = reshape(Jact,  [Njact_ex*Lact, Tall, Ntry]);