トポグラフィの作成

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;