DFSファイル
DFSファイルから三角メッシュと頂点を読み込む関数Mファイル
DFS file format
DFSファイル形式は以下のようになっています。
|
ヘッダー |
|
メッシュ数、頂点数、ポインタ |
|
メタデータ |
埋め込みXML |
不使用 |
|
被検者メタデータ |
埋め込みXML |
不使用 |
|
三角データ |
int32,int32,int32 |
メッシュ数×3×4バイト |
|
頂点座標 |
single,single,single |
0〜頂点数-1 |
|
頂点の法線 |
single,single,single |
オプション ないときは0 |
|
頂点のUV座標 |
single,single |
オプション テキストマップ座標 |
|
頂点のRGB色 |
single,single,single |
オプション 0〜1 |
|
頂点のラベル |
int16 |
オプション |
|
頂点属性 |
single |
オプション 曲率など |
とりあえず、MATLABで読んでみます。
'D:\MATLAB\akira\BrainSuite\akira.brain.dfs'
fid=fopen(loadname,'rb');
B=fread(fid,'*uint8');
fclose(fid);
ヘッダーの大きさとかどうなってるかよくわかりません。
B(1:32)'
char(B(1:32)')
結果です。
ans =
1×32 の uint8 行ベクトル
1 列から 16 列
68 70 83 95 76 69 32 118 50 46 48 0 184 0 0 0
17 列から 32 列
0 0 0 0 0 0 0 0 244 181 2 0 252 90 1 0
ans =
'DFS_LE v2.0 ¸
ôµ[1]
üZ '
>>
いろいろなdfsファイルみたんですが、どうも0〜23番地は同じ数字のようです。読み飛ばすことにします。
BrainSuiteによると
akira.brain.dfs
D:/MATLAB/akira/BrainSuite
triangles: 177652
vertices: 88828
normals: 88828
uv map: 0
vertex color: 0
vertex labels: 0
です。
24〜55番地を見ていきます。
address=(1:32)+24;
B(address)'
結果です。
ans =
1×32 の uint8 行ベクトル
1 列から 16 列
244 181 2 0 252 90 1 0 0 0 0 0 0 0 0 0
17 列から 32 列
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
>>
三角メッシュ数、頂点数だけ表示されていて、ポインタ部分は全部ゼロになっています。
とりあえず、以下を実行します。
typecast(B(25:28),'int32')
typecast(B(29:32),'int32')
結果です。
ans =
int32
177652
ans =
int32
88828
>>
メッシュ数と頂点数はBrainSuiteの情報と一致しています。
24番地がメッシュ数NT、28番地が頂点数ということでいいようです。
32番地から223番地を見ていきます。
address=(1:192)+32;
B(address)'
結果です。
ans =
1×192 の uint8 行ベクトル
1 列から 16 列
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
17 列から 32 列
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
33 列から 48 列
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
49 列から 64 列
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
65 列から 80 列
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
81 列から 96 列
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
97 列から 112 列
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
113 列から 128 列
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
129 列から 144 列
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
145 列から 160 列
0 0 0 0 0 0 0 0 0 0 0 0 236 1 0 0
161 列から 176 列
229 1 0 0 0 0 0 0 230 1 0 0 1 0 0 0
177 列から 192 列
229 1 0 0 230 1 0 0 0 0 0 0 1 0 0 0
>>
189番地から数字が出てきました。
ヘッダー部分の「,」にみえた184というのはヘッダー部分の大きさ、三角メッシュ情報のポインタだと勝手に解釈します。
以後0番目の頂点からなる三角メッシュのリストとしてデータを見ていきます。
nt=typecast(B(25:28),'int32');
nv=typecast(B(29:32),'int32');
address=184+(1:(3*nt*4));
tri=typecast(B(address),'int32');
tri=reshape(tri,[3,nt]);
max(tri(:))
結果です。
ans =
int32
88827
>>
頂点の番号は0〜88827までですのでこれでOKです。
次は頂点座標です。MATLABでは頂点0番は認識しないので1ずつ値を増やし、頂点番号を1からとします。
address=184+3*nt*4+(1:(3*nv*4));
vertex=typecast(B(address),'single');
vertex=reshape(vertex,[3,nv]);
figure;
faces=tri+1;
hbrain=patch('faces',faces','vertices',vertex');
set(hbrain,'facecolor',[1,1,0],'edgecolor','none');
axis tight;
daspect([1,1,1]);
view([-120,20]);
light('position',[-1,0,0]);
rotate3d on;
結果です。
BrainSuiteの情報では法線ベクトルもあるそうです。最後の行はホントに法線なのか１になるかどうかの確認です。
address=184+3*nt*4+3*nv*4+(1:(3*nv*4));
vnorm=typecast(B(address),'single');
vnorm=reshape(vnorm,[3,nv]);
max(sum(vnorm.^2,1))
ん？エラーが出ました。
インデックスが配列要素数 (3197944) を超えています。
エラー: Untitled (line 37)
vnorm=typecast(B(address),'single');
>>
今まで読み込んだファイル上のアドレスとファイルのデータサイズの長さとを比較してみます
184+3*nt*4+3*nv*4,length(B)
結果です。
ans =
int32
3197944
ans =
3197944
>>
頂点座標を読み込んだところで丁度ファイルの長さが終わっていました。
DFSファイルから三角メッシュと頂点を読み込む関数Mファイル
以下のような関数Mファイルを作成しました。
function R=hns_loaddfs(loadname)
R=[];
try
fid=fopen(loadname,'rb');
B=fread(fid,'*uint8');
fclose(fid);
catch
return
end
try
ptrimesh=typecast(B(13:16),'int32');
ntri=typecast(B(25:28),'int32');
nvertex=typecast(B(29:32),'int32');
pvertex=ptrimesh+3*ntri*4;
pnormal=typecast(B(33:36),'int32');
puv=typecast(B(37:40),'int32');
prgb=typecast(B(41:44),'int32');
plabel=typecast(B(45:48),'int32');
pattr=typecast(B(49:52),'int32');
catch
return
end
try
address=ptrimesh+(1:(3*ntri*4));
x=typecast(B(address),'int32');
x=reshape(x,[3,ntri]);
R.faces=x+1;
address=pvertex+(1:(3*nvertex*4));
x=typecast(B(address),'single');
R.vertices=reshape(x,[3,nvertex]);
catch
return
end
if pnormal>0
try
address=pnormal+(1:(3*nvertex*4));
x=typecast(B(address),'single');
R.normals=reshape(x,[3,nvertex]);
catch
return
end
end
if puv>0
try
address=puv+(1:(2*nvertex*4));
x=typecast(B(address),'single');
R.uv=reshape(x,[2,nvertex]);
catch
return
end
end
if prgb>0
try
address=prgb+(1:(3*nvertex*4));
x=typecast(B(address),'single');
R.rgb=reshape(x,[3,nvertex]);
catch
return
end
end
if plabel>0
try
address=plabel+(1:(nvertex*2));
x=typecast(B(address),'int16');
R.label=reshape(x,[3,nvertex]);
catch
return
end
end
if pattr>0
try
address=pattr+(1:(nvertex*4));
R.attr=typecast(B(address),'single');
catch
return
end
end
以下を実行しました。left/right.pial.cortex.dfsは色情報RGBがありましたので、その情報を使っています。
DFSfiles={
'D:\MATLAB\akira\BrainSuite\akira.left.pial.cortex.dfs'
'D:\MATLAB\akira\BrainSuite\akira.right.pial.cortex.dfs'
'D:\MATLAB\akira\BrainSuite\akira.scalp.dfs'
};
for n=1:length(DFSfiles)
V=hns_loaddfs(DFSfiles{n});
if n==1
figure('color',[1,1,1]);
end
h=patch('faces',V.faces','vertices',V.vertices',...
'FaceColor',[1,0.5,0],'FaceAlpha',0.5,'edgecolor','none');
if isfield(V,'rgb')
set(h,'FaceVertexCData',V.rgb','FaceColor','interp','FaceAlpha',1);
end
if n==1
hold
on;
daspect([1,1,1]);
end
end
axis tight;
view([-120,20]);
rotate3d on;
結果です。
GIfTIファイル化
GIfTI toolbox for MATLABを使います。左半球を表示させ、BrainVISA風のLhemi.giiという名前で保存します。
close all;
addpath('D:\MATLAB\gifti-1.6');
loadname='D:\MATLAB\akira\BrainSuite\akira.left.pial.cortex.dfs';
v=hns_loaddfs(loadname);
v.faces=v.faces';
v.vertices=v.vertices';
g=gifti(v);
figure('color',[1,1,1]);
plot(g);
save(g,'D:\MATLAB\akira\BrainSuite\Lhemi.gii','Base64Binary');
結果です。
警告: cameramenu
は将来のリリースで削除される予定です。代わりに cameratoolbar
を使用してください。
>>
このGIfTIファイルはFreeSurferのmris_convertコマンドを使うことでFreeSurferのsurfファイルに変換できます。
メッシュを膨らませる
パット見、BrainSuiteにはinflated brainの機能は見当たりません。
BrainVISAのMorphologistのsurface→Inflated Cortical Surface機能を使ってみます。
以下を実行しました。BrainVISAに合わせ、Lwhite.giiという名前で保存します。
close all;
addpath('D:\MATLAB\gifti-1.6');
loadname='D:\MATLAB\akira\BrainSuite\akira.left.inner.cortex.dfs';
v=hns_loaddfs(loadname);
v.faces=v.faces';
v.vertices=v.vertices';
g=gifti(v);
figure('color',[1,1,1]);
plot(g);
save(g,'D:\MATLAB\akira\BrainSuite\Lwhite.gii','Base64Binary');
Ubuntu 12.04の/home/ns/brainsuiteにコピーします。
BrainVISA 4.6.1を起動し、Morphologist→surfacd→Inflate Cortical Surfaceを選択します。
inputファイルを選択します。
/home/ns/brainsuite/Lwhite.giiを選択します。
output_meshのファイル名をxx_inflated.giiとし、Runを押します。
curvature textureのところが空欄になっていると実行しないので手入力します。
結果です。
M
filename2='C:\akira\BrainSuite\akira_Lwhite_inflated.gii';
filename3='C:\akira\BrainSuite\akira_Lwhite_curv.gii';
f2=gifti(filename2);
f3=gifti(filename3);
figure;
h=trimesh(f2.faces,f2.vertices(:,1),f2.vertices(:,2),f2.vertices(:,3),f3.cdata);
set(h,'edgecolor','none','facecolor','interp');
map=-gray(2)*0.5+0.75;
colormap(map);
set(gca,'clim',[-1,1]);
daspect([1,1,1]);
view([-120,20]);
rotate3d on;
