=====  VBMEG リリースノート

バージョン：0.50a_0
CVSタグ：version_0_5_0_a_0

=====  CVS情報

CVSROOT
  :pserver:username@lobule2:/home/cnb/rhayashi/cvs_repository/vbmeg_genesis


====== 前バージョン(0.41)からの修正および変更

1.VB推定の変更
2.ソフト制約条件の追加
3.脳モデル作成に関する変更
4.脳モデルの座標変換の修正 
5.GUIの変更
6.job_current_trの仕様変更
7.その他変更

============ 1.VB推定の変更

■推定方法修正のまとめ (job_vb)  ver 0.5

１．ベースライン活動強度推定の変更
・bayes_parm.forward_model = 'focal' の時
    ベースライン活動強度を推定する時、global_region = focal_region とした。
    VB推定時にはglobal_regionは無いので、focal_regionを使う方が妥当

・bayes_parm.forward_model = 'focal+global' の時
    global_region から focal_region を除く
    global_region と focal_region の干渉を避ける

・bayes_parm.remove_crossed_area = ON で以上の変更を実行
  bayes_parm.remove_crossed_area = OFF で、従来動作

・以上の変更をするために、vb_leadfield_preparationを
  vb_focal_leadfield_preparationとvb_global_leadfield_preparationに分けた

・以上の変更に対応して、vbmeg_current_reconstruct も修正

２．事前分散パラメタの自動推定
・bayes_parm.a0_act = [] の時、
  a0_actを Bayes Wiener 推定(vbmeg_multi_fmri)で自動推定するようにした

・a0_act に値が指定された時は、従来と同じ

・第1ステップ：
  fMRIパターンの係数を 各時間窓毎に Bayes Wiener 推定で求める

・第2ステップ：
  第1ステップで計算した電流分散を事前情報として用いてVB推定

・複数のfMRIパターンを指定可能
  複数のfMRIパターンの重ね合わせで電流分散を表しその重みを推定
  alpha = v0 + sum_n a0(n) * fMRI(n)^2
  として、ベースライン活動 v0 は固定して、a0(n) を推定
  各時間窓毎でこの重みは異なる

３．ソフト制約条件を使うVBMEG推定を実装
  詳細は別に記す

・内部関数 cortex_normal_variance でBVモデルでの近傍法線方向分散計算

４．VB推定の更新則の修正

・最初にVB更新則を使い、後半に高速更新則を使うように変更

bayes_parm.Npre_train : 最初のVB更新則を使う回数
bayes_parm.Ntrain     : 総更新回数

■各関数の主な変更点

job_vb
% Wiener estimation for hyper-prior parameter 'a0_act'
% Set  'global area' = 'focal area' in 'focal' model
% Remove focal area from global area

vb_set_vb_parm
% VB estimation parameteres 'vb_parm' are set based on 'bayes_parm'
% Input arguments are modified

set_multiple_prior
% Set (multiple) fMRI activity pattern(s) for Bayes Wiener estimate

set_vb_prior
% Set prior hyperparameter using fMRI activity pattern  or Wiener estimate

vbmeg_multi_fmri
% Bayesian Wiener estimation for MEG. 
% Current variance is represented by mixture of fMRI patterns
% Weight coefficient for each fMRI pattern is estimated

vbmeg_ard_estimate2
% VB-update equation is used for iteration <= Npre_train
% Soft normal constraint is implemented

vb_baseline_variance_acquisition
% vbmeg_multi_fmri is used for baseline estimation


============ 2.ソフト制約条件の追加

■ソフト制約条件の使い方 ver 0.5 ( For advanced user)

■ソフト制約条件用リードフィールドの作成

・法線平均計算フラグ average mode for normal direction
basis_parm.normal_mode
    = NORMAL_ADJACENT_MEAN = 0 : 隣接点の法線平均(従来計算法)
    = NORMAL_NEIGHBOR_MEAN = 1 : BV脳モデル近傍点の法線平均

・ソフト制約条件でVB推定するためには、以下条件でリードフィールドを作成
・basis_parm.Basis_mode  = CURRENT_3D % 独立な法線方向は3方向
・basis_parm.normal_mode = NORMAL_NEIGHBOR_MEAN

■ソフト制約条件を使うVBMEG推定

・ソフト制約条件フラグ : bayes_parm.soft_mode

・ハード制約条件 (従来法 : Default)
bayes_parm.soft_mode = 0;  
bayes_parm.variance_orientation = OFF;  

・ソフト制約条件指定
bayes_parm.soft_mode = 1;
bayes_parm.variance_orientation = ON;  

・ソフト制約条件を使う時は原則として低解像度脳モデルを使用する
  ただし、プログラム上は標準の脳モデルを使っても問題ない

■ソフト制約条件のリードフィールド整合性チェック

・bayes_parm.soft_mode = 1 の時、GUI で整合性をチェック
basis_parm.Basis_mode = 3 ＆ basis_parm.normal_mode = 1 を チェック
[Norient, normal_mode] = load_basis_info(basisfile);

■ソフト制約条件の方向分散パラメタ
結果が思わしくない時はこれを変えて試す

% tangential variance coefficient
bayes_parm.tan_var = 1 (Default)
% tangential variance max
bayes_parm.var_max = 1/sqrt(2) (Default)
% 法線方向分散 を 1 とした時、
% 接線方向分散 = sin( tan_var * 平均法線方向からの角度平均) <= var_max


============ 3.脳モデル作成に関する変更

■ job_brain, job_leadfield の修正点  ver 0.5

１．job_brainでの処理の追加

・ブレインボイジャー(Freesurfer)脳モデルの各頂点に皮質上で最も近い、
  脳モデル座標点を求め、brain_fileに変数追加

normal_stat = original_normal_statics(proj_root,brain_parm);

% normal_stat.neighbor_org : original BV brain neighbor index list
% normal_stat.normal_org   : original BV brain normal vector

・解像度を１／３にした低解像度脳モデルを作成
  最近傍点間の平均距離が約5mm 、頂点数約8000
  job_reduced_cortex(proj_root,brain_parm)

・１／３低解像度脳モデルのファイル名
  元の brain_file  = [ xxx.brain.mat ]
  低解像度モデル   = [ xxx_reduce.brain.mat ]
  area_file, act_file も同様

・job_inflate
brain_file に対するinflate モデルの作成と
低解像度脳モデルに対するinflate モデルの作成を行う

・追加パラメタ

% BV脳モデルの各点に最も近い皮質頂点を探す時の最大半径
% Max radius for BV surface neighbor
brain_parm.R_normal = 0.004; % Default 4mm
% Reduce ratio for reduced brain
brain_parm.reduce_brain = 1/3;  % Default = 1/3;

２．リードフィールド計算方法の追加

・MEGリードフィールドを球面調和関数展開で計算
・EEGリードフィールドを球モデルで計算
・EEGリードフィールドをBEMで計算

・追加パラメタ

basis_parmに以下のパラメタを追加したが、通常の場合これらは設定しなくても
プログラム中でデフォルト設定を使うので問題ない。


・BEMやMEG球面調和関数計算の時
% basis_parm.Recalc     : Recalculation force flag    (Default OFF)
    = ON の時は、ファイルに MEGinfo.SPHinfo や Omega が存在しても
      これらを再計算する

・MEG球面調和関数計算の時
% basis_parm.func_order : Order of spherical harmonics function (Default 35)

・EEG
% basis_parm.sigma      : Conductivity of each region
% basis_parm.radius     : Relative radius of sphere
basis_parm.sigma = [ 1.0  0.0125  1.0 ]; (Default)
basis_parm.radius = [ 0.87  0.92   1.0 ];(Default)

３．既存の脳モデルへの法線方向情報の追加

・既存の brainfile に法線方向情報を追加し、低解像度脳モデルを作る関数

job_normal_stat(proj_root,brain_parm)

４．領域情報(areakey)、活動強度情報(actkey) を低解像度脳モデル用に変換する

・ 領域情報(areakey)、活動強度情報(actkey) を低解像度脳モデル用に変換する関数
% acttivity map is transformed to act map in the reduced brain model
act_to_reduce_brain(actfile, areafile, actkey)
act_to_reduce_brain(actfile, areafile, actkey, reducekey,reduced_actfile)

% Area index is transformed to area index in the reduced brain model
area_to_reduce_brain(areafile, areakey)
area_to_reduce_brain(areafile, areakey, reducekey,reduced_areafile)


============ 4.脳モデル作成に関する変更

■脳モデルの座標変換の修正 ver 0.5

これまでBrainVoyagerからSPM座標への変換は島津の梶原さんが作ったコードを
もとに行っていたが、以下のようにボクセル座標とmm実座標との変換で不整合が
見つかったので、修正した。

このため、ver 0.4 以前に計算した脳モデルと ver 0.5 以降に計算した脳モデルの
座標は少し異なっている。

■データの整合性に関する注意点

１．ver 0.4 以前に計算した脳モデルを ver 0.5 で使う場合
脳モデル1式(*.brain.mat, *.act.mat, *.area.mat, *.basis.mat, *.meg.mat)
を持っている場合は、これをそのまま使ってjob_vb, job_current を行える。

２．ver 0.5 で最初から脳モデルを作り直して解析する(推奨)

３．注意点
ver 0.4 で作った*.brain.matに対応する*.act.mat, *.area.mat と
ver 0.5 で作った*.brain.matに対応する*.act.mat, *.area.mat は異なるので、
*.brain.matを作り直したら、*.act.mat, *.area.matも作り直さなければいけない。

% ---- make_brain_data2 (ver 0.5) ---
Brain Voyager coordinate(V) から Right-hand SPM coordinate(Vspm) へ
直接変換する関数 bvoyger_spm_right を作った。
V の原点[255/2. 255/2. 255/2] と Vspm の原点[ 0, 0, 0]を合わせる

  Vspm = bvoyger_spm_right(V,2)
    Vspm(:,1) = 255/2 - V(:,3); 
    Vspm(:,2) = 255/2 - V(:,1); 
    Vspm(:,3) = 255/2 - V(:,2); 

% ---- make_brain_data ( before ver 0.5) ---
これまでのプログラムでは、以下のように2段階の座標変換を行っていた。

Vox  = bvoyger_to_vox(V)    
Vspm = vox_to_spm_right(Vox)

これは以下の座標変換に対応する
    Vspm(:,1) = 255/2 - V(:,3) + 1  ;
    Vspm(:,2) = 255/2 - V(:,1) - 1/2;
    Vspm(:,3) = 255/2 - V(:,2) - 1/2;

すなわち、直接原点同士を合わせるのと0.5-1.0mm ずれている。

<参考>
Vox = bvoyger_to_vox(V)
    Vox(:,1) = V(:,1) + 1; 
    Vox(:,2) = V(:,2) + 1; 
    Vox(:,3) = 256 - V(:,3) - fix( (256 - Vdim(1)*Vsize(1))/2 );

Vspm = vox_to_spm_right(Vox)
    Vspm(:,1) =   Vox(:,3) - Vdim(1)*Vsize(1)*0.5;
              = - V(:,3)   + 255/2 + 1
    Vspm(:,2) = - Vox(:,1) + Vdim(2)*Vsize(2)*0.5;
              = - V(:,1)   + Vdim(2)*Vsize(2)/2 
              = - V(:,1)   + 255/2 - 1/2
    Vspm(:,3) = - Vox(:,2) + Vdim(3)*Vsize(3)*0.5;
              = - V(:,2)   + Vdim(3)*Vsize(3)/2
              = - V(:,2)   + 255/2 - 1/2

============ 5.GUIの変更

1.GUIのリソースをスクリプト管理していたものをfigureファイル管理に変更した。
  (前処理GUIのAdvancedパラメータGUIは対応未完了。今後置き換える)
2.推定GUIの変更点
  1.以下のパラメータ追加に対応した(各パラメータについては、「1.VB推定の変更」を参照)
     bayes_parm.Npre_train
               .soft_mode
               .tan_var
               .var_max
               .removed_cross_area
  2.以下のパラメータ変更に対応した
     bayes_parm.a0_act(パラメータ横にautoのチェックボックスを表示し、
                       autoにチェックが入っているときは、a0_act=[]としてjob_vbを呼ぶ)

  3.デフォルト値の取得を, template/set_bayes_parm.mから取得するよう変更
  4.仕様変更
    (a)推定モデルをFoward modelとnoise estimation modelからなる構成に変更
    (b)act keyを複数選択できるよう仕様変更
  5.ノイズ推定モデルのTimeWindowをファイルから読み込むことができない不具合を修正した

3.前処理GUIの変更点
  - Leadfield計算GUI
    ・Bem modeにBASIS_MEG_HARMONICSを追加した
    ・Normal modeを追加した
4.project_mgrの履歴パラメータの再利用機能追加
  - project_mgrに表示されている、実行パラメータ履歴の１つをクリックし、
    右クリックを押下すると、メニューが表示され、launch applicationを選択すると、
    選択したパラメータを持った状態でGUI(前処理・推定とも)が起動する機能を追加した。

============ 6.job_current_trの仕様変更

■  job_current_trの仕様変更

・常にトライアル毎の電流を計算するように変更
  parm.trial_average = OFF
  トライアル平均電流のみを計算する時は job_current を使う。
・jactdir の下にトライアル毎の電流ファイルが保存される
・currfile にはトライアル平均電流とトライアル毎の電流に関する情報が保存される
・job_current を trial_average = OFF で実行すると
  トライアル毎の電流がひとつの currfile に保存される

job_current_tr(proj_root,parm)

--- Input parameters
 - parm.resultfile : Estimated result by job_vb
   .bayes.mat file.
 - parm.currfile   : Average current is saved.
                     Infomation of current file are also saved.
                     current of each trials are saved into 'jactdir' directory
   .curr.mat file 
 - parm.jactdir : output directory of current files for each trial
                  relative path from proj_root
                  do not include '/' or '\' in the end of string
                  If 'jactdir' is not given,  jactdir = 'temp'
--- Optional  (default)
 - parm.areafile ([])  --> focal window
   .area.mat file.
 - parm.area_key ([])
   ID of interested area.
 - parm.tsubsmpl ([]) --> whole time points
   Time points where currents are calculated. For example,
   if tsubsmpl = [1:5:100], currents are calculated at only every 5 sample
   point.
--- Variables saved in output file (.curr.mat)

 - Jact
 - Jbck
 - ix_act
 - ix_bck
 - Lact
 - Lbck
 - Varact
 - bayes_parm
 - vb_parm

============ 7.その他変更

■バグ修正
1.job/job_inflate.m
  ディレクトリパスの設定に不具合があり、正常なパラメータで実行しても
  check_inflate_modelでエラー終了してしまうバグを修正。
2.job/job_preprocess_meg3.m
  内部関数save_divisionにあったバグを修正
3.job/job_current_tr.m, job_current.m
  - 入力パラメータのチェックに関するバグ修正
  - Tstartに関するバグ修正

■関数追加

1.vbmeg/common/utility/kl_distance.m
  kl_disntace関数追加

  KL-distance of X and Y
  KL = kl_distance(X,Y)
  X and Y are normalized before calculation 
  Zero value component is replaced by eps

2.common/loadfunc/load_basis_info.m
  ・Basisファイル情報読み込み関数追加

  [Norient, normal_mode] = load_basis_info(basisfile);
  [Norient, normal_mode, bem_mode] = load_basis_info(basisfile);
  --- Input
  basisfile : leadfileld basis file name
  --- Output
  Norient     : number of independent current orientation 
  normal_mode : normal average mode
  basis_parm  : leadfield calculation parameters

3.common/loadfunc/is_inflate_cortex.m
  インフレートデータがbrainファイルに存在するかチェックする関数を追加
  result = is_inflate_cortex(brainfile)
 
  -- Input 
  brainfile   : Brain file name
  --- Output
    result = SUCCESS;
    result = FAILURE;

4.元のBrainVoyagerの法線情報があるかどうかチェックする関数を追加
  common/loadfunc/is_normal_statics.m 

  Check existence of original BV brain normal information
   result = is_normal_statics(brainfile)
  --- Output
    result = SUCCESS;
    result = FAILURE;

5.Atlas labelを取得するための関数を追加
  tool_box/atlas2vb/get_atlas_area_label.m

  Get area index corresponding to atlas map in the given area
 
  [ix_list, ix_name] = get_atlas_label(ix_area,atlasfile,atlas_id)
  --- Input
  ix_area : Vertex index for given area
  atlasfile : File name of subject atlas file (*.atlas.mat)
  atlas_id  : Atlas name ('aal', 'brodmann')
  --- Output
  ix_list : set of area index corresponding to atlas region
  ix_name : name of atlas region corresponding to 'ix_list'

6.領域情報(areakey)、活動強度情報(actkey) を低解像度脳モデル用に変換する関数
  % acttivity map is transformed to act map in the reduced brain model
  act_to_reduce_brain(actfile, areafile, actkey)
  act_to_reduce_brain(actfile, areafile, actkey, reducekey,reduced_actfile)

  % Area index is transformed to area index in the reduced brain model
  area_to_reduce_brain(areafile, areakey)
  area_to_reduce_brain(areafile, areakey, reducekey,reduced_areafile)

7.既存の brainfile に法線方向情報を追加し、低解像度脳モデルを作る関数
  job/job_normal_stat.m

  job_normal_stat(proj_root,brain_parm)

8.common/utility/matlab_version.m
  MATLABのバージョン取得を行う関数を追加。

  Usage: [major, minor, revision, build] = matlab_version;
  [OUT]
      major : MATLAB major version
      minor : MATLAB minor version
      revision : MATLAB revision
      build : MATLAB build number

9.template/
  以下のファイルを追加。
  set_current_parm_test.m, set_current_parm.m (job_current用パラメータファイル)
  set_bayes_default_parameters.m

■変更・修正
1.MATLABのバージョンをチェックして実行している箇所の修正
  ※MATLABのバージョン取得をMATLAB組み込み関数version('-release')から
    関数：matlab_versionを使用するように変更。
    MATLABが返却するリリース番号の仕様が大きく変更されたため。
  修正ファイル：vb_save, vb_fsave, dilation_3d.m erosion_3d.m test_morfology.m
                create_gui.m, proj_run.m, 
2.common/utility/err_array.m
  どんな種類の配列でも比較できるようにした
3.common/container/add_area.m, add_act.m
  データの保存をsave関数ではなく、vb_fsave関数を使用するよう変更.
4.plotfunc/check_model/
  関数コール時に必要なデータが入力されなかった場合は、以下の
  メッセージを表示してリターンするように変更。
    check_sphere_slice.m
        There is no analize file: *****
        Check of sphere model is skipped
    check_brain_modelm
        There is no analize file: *****
        Check of MRI and brain model is skipped
    check_fmri.m
        There is no inflate data: *****
        Plot of fMRI on inflate brain model is skipped
5.common/loadfunc/load_normal_statics.m
  brainファイル中に'normal_stat'が無い状態で、関数が呼ばれたときのエラーチェックを追加
6.vbmeg.mを修正
  ・旧バージョンのGUIに関する記述を削除
  ・バージョン表示機能を追加。書式：vbmeg('-version');

