find minimum path along cortical surface
[ID, DD] = vb_find_geodesic(Jstart, Jend, xxF, xxD ,Rmax)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
最短距離経路探索
隣接近傍点リスト( xxF ) とその距離( xxD ) を使って
皮質に沿って最短距離で開始点から終点までの経路を探索
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
INPUT
Jstart : 開始点のインデックス
Jend : 終点のインデックス
xxF{i} : 隣接近傍点インデックス
xxD{i} : 隣接近傍点距離
Rmax : 探索最大半径
OUTPUT
ID : 経路点インデックス
DD : 累積距離
Made by Masa-aki Sato on 2005-1-28
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function [ID, DD] = vb_find_geodesic(Jstart, Jend, xxF, xxD ,Rmax) 0002 % find minimum path along cortical surface 0003 % [ID, DD] = vb_find_geodesic(Jstart, Jend, xxF, xxD ,Rmax) 0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0005 % 最短距離経路探索 0006 % 隣接近傍点リスト( xxF ) とその距離( xxD ) を使って 0007 % 皮質に沿って最短距離で開始点から終点までの経路を探索 0008 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0009 % INPUT 0010 % Jstart : 開始点のインデックス 0011 % Jend : 終点のインデックス 0012 % xxF{i} : 隣接近傍点インデックス 0013 % xxD{i} : 隣接近傍点距離 0014 % Rmax : 探索最大半径 0015 % OUTPUT 0016 % ID : 経路点インデックス 0017 % DD : 累積距離 0018 % 0019 % Made by Masa-aki Sato on 2005-1-28 0020 % 0021 % Copyright (C) 2011, ATR All Rights Reserved. 0022 % License : New BSD License(see VBMEG_LICENSE.txt) 0023 0024 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0025 % アルゴリズム ： 最近傍点リストを用いた Tree 探索 0026 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0027 %１．自分自身をルートにする 0028 %２．各ルート点の隣接近傍点の集合を候補点リストに登録する 0029 %３．候補点リストの中でこれまでに試されていない新しい候補点を抽出する 0030 %４．出発点から候補点までの累積距離を計算(経路に依存) 0031 %５．出発点から候補点までの最短累積距離を求める 0032 %６．累積距離が指定半径以内にある点を見つける 0033 %７．上記の点を新しいルートにして２．へ戻る 0034 %８．終点に到達したら終了 0035 %９．指定半径以内にある新しいルートが無くなったときに終了 0036 0037 nextID = Jstart; % 近傍点インデックスリスト 0038 prevID = 0; % 1つ前の点のインデックスリスト 0039 nextDD = 0; % 近傍点累積距離リスト 0040 rootix = Jstart; % ルートのインデックスリスト 0041 rootd = 0; % ルートの累積距離リスト 0042 itree = 0; % tree number 0043 0044 while itree==0 | ~isempty(rootix), 0045 itree = itree+1; 0046 nroot = size(rootix,1); % 今回のルート点数 0047 0048 % 候補点 : 全てのルートの隣接近傍点 0049 ixlist = []; % 隣接近傍点インデックスリスト 0050 jplist = []; % 1つ前の点のインデックス 0051 ddlist = []; % 隣接近傍点への累積距離リスト 0052 0053 for i=1:nroot, 0054 iroot = rootix(i); % ルートのインデックス 0055 ixlist = [ixlist; xxF{iroot}]; % 隣接近傍点インデックス 0056 ddlist = [ddlist; xxD{iroot} + rootd(i) ]; % 隣接近傍点への累積距離 0057 newnum = length(xxF{iroot}); 0058 jplist = [jplist; iroot(ones(newnum,1))]; % 1つ前の点のインデックス 0059 end; 0060 0061 % 候補点の絞り込み 0062 ixuniq = unique(ixlist); % 重複インデックスを削除 0063 nextix = vb_setdiff2(ixuniq,nextID); % 探索済みのインデックスを削除 0064 0065 % 探索候補点インデックス 0066 nextix = nextix(:); 0067 Nnext = length(nextix); 0068 nextd = zeros(Nnext,1); 0069 prevj = zeros(Nnext,1); 0070 0071 for i=1:Nnext, % Loop of new dipoles 0072 % ixlist の中で nextix(i) に等しいインデックスを抽出 0073 % Jstart から nextix(i) へ複数の経路に対応 0074 jx = find( nextix(i)==ixlist ); 0075 0076 % nextix(i) への最短累積距離を探す 0077 [dmin ,jmin] = min(ddlist(jx)); 0078 % 最短累積距離 0079 nextd(i) = dmin; 0080 % 1つ前の点のインデックス 0081 prevj(i) = jplist(jx(jmin)); 0082 end 0083 0084 % 終点に到達したかチェック 0085 jfinal = find(nextix == Jend); 0086 0087 if isempty(jfinal) 0088 % 最大半径以内の候補点を探す 0089 okix = find(nextd<Rmax); % 最大半径以内の候補点を探す 0090 rootix = nextix(okix); % 最大半径以内の候補点(次回のルート点) 0091 rootd = nextd(okix); % その最短累積距離 0092 prevj = prevj(okix); % 1つ前の点のインデックス 0093 0094 % 近傍点リスト 0095 nextID = [nextID ; rootix ];% 近傍点リスト 0096 prevID = [prevID ; prevj ]; % 1つ前の点のインデックス 0097 nextDD = [nextDD ; rootd ];% 最短累積距離 0098 else 0099 nextID = [nextID ; Jend ];% 近傍点リスト 0100 prevID = [prevID ; prevj(jfinal)]; % 1つ前の点のインデックス 0101 nextDD = [nextDD ; nextd(jfinal)]; % 最短累積距離 0102 break; 0103 end 0104 end 0105 0106 ID = []; 0107 DD = []; 0108 0109 Nend = length(nextID); 0110 iend = Nend; 0111 % Backtrack 0112 % 1つ前の点のインデックスをもとに終点から開始点へ経路探索 0113 0114 for n=1:Nend 0115 ID = [ID; nextID(iend)]; 0116 DD = [DD; nextDD(iend)]; 0117 jpre = prevID(iend); 0118 0119 if jpre == 0, break; end; 0120 0121 iend = find( jpre == nextID ); 0122 end 0123 0124 ID = flipud(ID); 0125 DD = flipud(DD);