トポグラフィの作成
MATLAB®には簡単にトポグラフィを作るための関数が充実しています。
平成24年9月20日の広大霞キャンパス基礎棟2Fで行ったMATLABセミナーの資料の一部です。 参加人数は7人くらいでした。
Contents
国際10-20法の電極位置
脳波の10-20法の電極位置座標を決定します。極座標を使います。F3,F4,P3,P4は適当です。
EEG=zeros(21,3);% 21個の電極のx,y,z座標 for n=1:5 switch n case 1;theta=0;phi=pi/2; ch=1;% Cz case 2;theta=(0:0.5:1.5)*pi;phi=0.6*pi/2; ch=2:5;% C4,Fz,C3,Pz case 3;theta=(0.25:0.5:1.75)*pi;phi=0.4*pi/2; ch=6:9;% F4,F3,P3,P4 case 4;theta=(0:0.2:1.8)*pi;phi=0.2*pi/2; ch=10:19;% T4 F8 Fp2 Fp1 F7 T3 T5 O1 O2 T6 case 5;theta=(0:1)*pi;phi=0; ch=20:21;% A2 A1 end; nch=length(ch);r=cos(phi); EEG(ch,:)=[r*cos(theta);r*sin(theta);sin(phi)*ones(1,nch)]'; end;
国際10-20法電極名の並び替え
チャンネル順はCz,C4,Fz,・・・です。 たぶん実際の脳波データはFp1,Fp2,F3、・・・の順になってます。 strcmpi関数を使って脳波データのチャンネル順になるようにの順番を並び替えます。
ChName={'Cz','C4','Fz','C3','Pz',... 'F4','F3','P3','P4','T4',... 'F8','Fp2','Fp1','F7','T3',... 'T5','O1','O2','T6','A2',... 'A1'};% 上で作った電極の順番 ChOrder={'Fp1','Fp2','F3','F4','C3','C4','P3','P4','O1','O2',... 'F7','F8','T3','T4','T5','T6','Fz','Cz','Pz','A1','A2'}; ChNo=zeros(length(ChOrder)); for n=1:length(ChOrder); ChNo(n,:)=strcmpi(ChName,ChOrder(n));% 大文字と小文字で区別しない end; ChNo=ChNo*diag(1:length(ChOrder)); ChNo=ChNo'; ChNo(ChNo==0)=[]; newEEG=EEG(ChNo,:);
電極位置表示
電極の位置を二次元と三次元で表示します。 三次元はマウスのドラッグで視点を変えることができます。
nch=size(newEEG,1); figure('color',[1,1,1]); subplot(1,2,1);% 二次元表示 theta=atan2(newEEG(:,2),newEEG(:,1)); phi=1-asin(newEEG(:,3))/pi*2; P=[cos(theta).*phi,sin(theta).*phi]; plot(P(:,1),P(:,2),'ro'); for n=1:nch; text(P(n,1),P(n,2),ChOrder(n)); end; hold on; for n=1:3 switch n case 1;r=1; case 2;r=0.8; case 3;r=0.4; end; h=rectangle('position',r*[-1,-1,2,2],'curvature',[1,1]); if n>1;set(h,'linestyle',':');end; end; patch([-0.1,0,0.1],[0.9,1,0.9]+0.1,[1,1,1]); axis off;title('二次元表示'); subplot(1,2,2);% 三次元表示 plot3(newEEG(:,1),newEEG(:,2),newEEG(:,3),'ro'); for n=1:nch; text(newEEG(n,1),newEEG(n,2),newEEG(n,3),ChOrder(n)); end; hold on; [x,y,z]=sphere(36);r=0.95; x=x(18:end,:)*r;y=y(18:end,:)*r;z=z(18:end,:)*r;%球の下半分は無視 h=surf(x,y,z); set(h,'facecolor',[1,1,1]*0.9,'linestyle','none'); n=[6,13];plot3(x(n,:)',y(n,:)',z(n,:)','k:'); axis off;title('三次元表示'); daspect([1,1,1]); rotate3d on;
mambrane関数の結果を脳波電極の値にして等電位地図の作成
ある一時点の各脳波電極の電位分布がmembrane関数になったと仮定し、 griddata関数を使って等電位地図を作成します。 具体的にはmembrane関数で得られた電位分布から各脳波電極の電位を取得し、 各電極の電位データからgriddata関数を使って電位分布を再合成して 等電位地図を作成します。contour関数は二次元、surfは三次元用です。
M=membrane(2);% membrane関数をトポグラフィにします % M=peaks(50);% peaks関数をトポグラフィにします Msize=size(M,1); [XI,YI]=meshgrid(1:Msize,1:Msize); newP=round((P+1)/2*(Msize-1))+1;% 電極位置をトポグラフィの座標に変換 newP(:,3)=0; nP=size(newP,1); for n=1:nP; newP(n,3)=M(newP(n,2),newP(n,1));% X-Y -> Y-X end; white=[1,1,1]; figure('color',white); subplot(2,2,1); contour(M);hold on; plot(newP(:,1),newP(:,2),'ko','MarkerFaceColor',white); rectangle('position',[0.5,0.5,Msize*[1,1]],'curvature',[1,1]); for n=1:nP;text(newP(n,1),newP(n,2),ChOrder(n));end; axis off; title('元の等電位図'); subplot(2,2,2); surf(M);hold on; plot3(newP(:,1),newP(:,2),newP(:,3),'ko','MarkerFaceColor',white); for n=1:nP;text(newP(n,1),newP(n,2),newP(n,3),ChOrder(n));end; axis tight;rotate3d on; title('元のトポグラフィ'); ZI=griddata(newP(:,1),newP(:,2),newP(:,3),XI,YI);% 強力トポグラフィ作成関数 subplot(2,2,3); contour(ZI);hold on; plot(newP(:,1),newP(:,2),'ko','MarkerFaceColor',white); rectangle('position',[0.5,0.5,Msize*[1,1]],'curvature',[1,1]); for n=1:nP;text(newP(n,1),newP(n,2),ChOrder(n));end; axis off; title('脳波の等電位図'); subplot(2,2,4); surf(ZI);hold on; plot3(newP(:,1),newP(:,2),newP(:,3),'ko','MarkerFaceColor',white); for n=1:nP;text(newP(n,1),newP(n,2),newP(n,3),ChOrder(n));end; axis tight;rotate3d on; title('脳波のトポグラフィ'); Lim=[0,Msize+1]; set(findobj(gcf,'type','axes'),'xlim',Lim,'ylim',Lim,'zlim',[-1,1]);
DELAUNAYで三次元トポグラフィに
delaunay関数を使うと三角メッシュが作成できます。 このメッシュとtrisurf関数を使うとopenGLを用いた綺麗な三次元トポグラフィが作成できます。
tri=delaunay(P(:,1),P(:,2));% 三角メッシュの作成 figure('color',[1,1,1]); h=trisurf(tri,newEEG(:,1),newEEG(:,2),newEEG(:,3),newP(:,3)); set(h,'facecolor','interp');hold on;axis off; for n=1:nch; text(newEEG(n,1),newEEG(n,2),newEEG(n,3),ChOrder(n)); end; daspect([1,1,1]);rotate3d on;
国際10-10法の電極位置
この部分はセミナーでは発表してません。追加したものです。 脳波の10-10法の電極位置座標を決定します。 電極の名称がT3→T7、T4→T8になります。またA1とA2は今回省略します。 最初に基準となる49個の電極位置を決定します。
EEG=zeros(49,2);% 個の電極のtheta,phi座標 for n=1:6 switch n case 1;theta=0;phi=pi/2; ch=1;% Cz case 2;theta=(0:0.5:1.5)*pi;phi=0.8*pi/2; ch=2:5;% C2,FCz,C1,CPz case 3;theta=(0:0.5:1.5)*pi;phi=0.6*pi/2; ch=6:9;% C4,Fz,C3,Pz case 4;theta=(0:0.5:1.5)*pi;phi=0.4*pi/2; ch=10:13;% C6,AFz,C5,POz case 5;theta=(0:0.1:1.9)*pi;phi=0.2*pi/2; ch=14:33;% T8,FT8,F8,AF8,Fp2,Fpz,Fp1,AF7,F7,FT7,T7,TP7,P7,PO7,O1,Oz,O2,PO8,P8,TP8 case 6;theta=[0,0.1,0.2,0.5,0.8:0.1:1.9]*pi;phi=0*pi/2; ch=34:49;% T10,FT10,F10,Nz,F9,FT9,T9,TP9,P9,PO9,l1,lz,l2,PO10,P10,TP10 end; nch=length(ch); EEG(ch,:)=[theta',ones(nch,1)*phi]; end; theta=EEG(:,1);phi=EEG(:,2);r=cos(phi); EEG=[r.*cos(theta),r.*sin(theta),sin(phi)]; % 次に残りの36個の電極中間点の電極位置を決定します。 EEG(50:85,:)=0; for n=1:7 switch n case 1;ch=50:53;% F4,F3,P3,P4 ch0=[16,7;7,22;9,26;32,9]; case 2;ch=54:61;% FC4,F2,F1,FC3,CP3,P1,P2,CP4 ch0=[6,50;50,7;7,51;51,8;8,52;52,9;9,53;53,6]; case 3;ch=62:69;% F6,AF4,AF3,F5,P5,P03,P04,P6 ch0=[16,50;18,50;20,51;11,51;26,52;28,52;30,53;32,53]; case 4;ch=70:73;% FC2 FC6 AF6 AF2 ch0=[3,54;54,15;17,63;63,11]; case 5;ch=74:77;% FC1,AF1,AF5,FC5 ch0=[57,3;11,64;64,21;23,57]; case 6;ch=78:81;% CP5,CP1,P05,PO1 ch0=[25,58;58,5;27,67;67,13]; case 7;ch=82:85;% PO2,P6,CP2,CP6 ch0=[13,68;68,31;5,61;61,33]; end; EEG(ch,:)=(EEG(ch0(:,1),:)+EEG(ch0(:,2),:))/2; end; % 50~85の電極を球面上になるよう補正します。 Ch=EEG(50:85,:); EEG(50:85,:)=Ch./(sum(Ch.^2,2).^0.5*[1,1,1]); newEEG=EEG; ChOrder={'Cz','C2','FCz','C1','CPz','C4','Fz','C3','Pz',... 'C6','AFz','C5','POz','T8','FT8','F8','AF8','Fp2',... 'Fpz','Fp1','AF7','F7','FT7','T7','TP7','P7','PO7',... 'O1','Oz','O2','PO8','P8','TP8','T10','FT10','F10',... 'Nz','F9','FT9','T9','TP9','P9','PO9','l1','lz','l2',... 'PO10','P10','TP10','F4','F3','P3','P4','FC4','F2','F1',... 'FC3','CP3','P1','P2','CP4','F6','AF4','AF3','F5','P5',... 'P03','P04','P6','FC2','FC6','AF6','AF2','FC1','AF1',... 'AF5','FC5','CP5','CP1','P05','PO1','CP5','CP1','P05','PO1'};
電極位置表示
nch=size(newEEG,1); figure('color',[1,1,1]); theta=atan2(newEEG(:,2),newEEG(:,1)); phi=1-asin(newEEG(:,3))/pi*2; P=[cos(theta).*phi,sin(theta).*phi]; plot(P(:,1),P(:,2),'ro'); for n=1:nch; text(P(n,1),P(n,2),ChOrder(n)); end; hold on; for n=1:3 switch n case 1;r=1; case 2;r=0.8; case 3;r=0.4; end; h=rectangle('position',r*[-1,-1,2,2],'curvature',[1,1]); if n>1;set(h,'linestyle',':');end; end; patch([-0.1,0,0.1],[0.9,1,0.9]+0.1,[1,1,1]); axis off;title('二次元表示'); figure('color',[1,1,1]); plot3(newEEG(:,1),newEEG(:,2),newEEG(:,3),'ro'); for n=1:nch; text(newEEG(n,1),newEEG(n,2),newEEG(n,3),ChOrder(n)); end; hold on; [x,y,z]=sphere(36);r=0.95; ps=16; x=x(ps:end,:)*r;y=y(ps:end,:)*r;z=z(ps:end,:)*r;%球の下半分は無視 h=surf(x,y,z); set(h,'facecolor',[1,1,1]*0.9,'linestyle','none'); n=[6,13];plot3(x(n,:)',y(n,:)',z(n,:)','k:'); axis off;title('三次元表示'); daspect([1,1,1]); rotate3d on;
mambrane関数の結果を脳波電極の値にして等電位地図の作成
M=membrane(2);% membrane関数をトポグラフィにします % M=peaks(50);% peaks関数をトポグラフィにします Msize=size(M,1); [XI,YI]=meshgrid(1:Msize,1:Msize); newP=round((P+1)/2*(Msize-1))+1;% 電極位置をトポグラフィの座標に変換 newP(:,3)=0; nP=size(newP,1); for n=1:nP; newP(n,3)=M(newP(n,2),newP(n,1));% X-Y -> Y-X end; white=[1,1,1]; figure('color',white); subplot(1,2,1); contour(M);hold on; plot(newP(:,1),newP(:,2),'ko','MarkerFaceColor',white); rectangle('position',[0.5,0.5,Msize*[1,1]],'curvature',[1,1]); for n=1:nP;text(newP(n,1),newP(n,2),ChOrder(n));end; axis off; ZI=griddata(newP(:,1),newP(:,2),newP(:,3),XI,YI);% 強力トポグラフィ作成関数 subplot(1,2,2); contour(ZI);hold on; plot(newP(:,1),newP(:,2),'ko','MarkerFaceColor',white); rectangle('position',[0.5,0.5,Msize*[1,1]],'curvature',[1,1]); for n=1:nP;text(newP(n,1),newP(n,2),ChOrder(n));end; axis off; title('脳波の等電位図'); figure('color',white); title('元の等電位図'); subplot(1,2,1); surf(M);hold on; plot3(newP(:,1),newP(:,2),newP(:,3),'ko','MarkerFaceColor',white); for n=1:nP;text(newP(n,1),newP(n,2),newP(n,3),ChOrder(n));end; axis tight;rotate3d on; title('元のトポグラフィ'); subplot(1,2,2); surf(ZI);hold on; plot3(newP(:,1),newP(:,2),newP(:,3),'ko','MarkerFaceColor',white); for n=1:nP;text(newP(n,1),newP(n,2),newP(n,3),ChOrder(n));end; axis tight;rotate3d on; title('脳波のトポグラフィ'); Lim=[0,Msize+1]; set(findobj(gcf,'type','axes'),'xlim',Lim,'ylim',Lim,'zlim',[-1,1]);
DELAUNAYで三次元トポグラフィに
tri=delaunay(P(:,1),P(:,2));% 三角メッシュの作成 figure('color',[1,1,1]); h=trisurf(tri,newEEG(:,1),newEEG(:,2),newEEG(:,3),newP(:,3)); set(h,'facecolor','interp');hold on;axis off; for n=1:nch; text(newEEG(n,1),newEEG(n,2),newEEG(n,3),ChOrder(n)); end; daspect([1,1,1]);rotate3d on;