FreeSurferのfile
format
何故かbig endianです。
mriのT1.mgz
以下を実行します。
mgh=gunzip(filename);
fid=fopen(mgh{1},'rb');
B=fread(fid,'*uint8');
fclose(fid);
delete(mgh{1});
F=[];
k=1:28;
x=swapbytes(typecast(B(k),'int32'));
F.version=x(1);
F.width=x(2);
F.height=x(3);
F.depth=x(4);
F.nframes=x(5);
F.type=x(6);% 0
uint8 1 int32 3 single 4 int16
switch F.type
case 0;fmt='uint8';fbyte=1;
case 1;fmt='int32';fbyte=4;
case 3;fmt='single';fbyte=4;
case 4;fmt='int16';fbyte=2;
end
F.dof=x(7);%
degree of freedom
k=29:30;
F.goodRASFlag=swapbytes(typecast(B(k),'int16'));
k=31:90;
x=swapbytes(typecast(B(k),'single'));
F.spacingx=x(1);
F.spacingy=x(2);
F.spacingz=x(3);
F.xr=x(4);%RAS
right anterior superior
F.xa=x(5);
F.xs=x(6);
F.yr=x(7);
F.ya=x(8);
F.ys=x(9);
F.zr=x(10);
F.za=x(11);
F.zs=x(12);
F.cr=x(13);
F.ca=x(14);
F.cs=x(15);
F
k=91:284;
B(k)'
k=284+(1:F.width*F.height*F.depth*fbyte);
F.data=swapbytes(typecast(B(k),fmt));
F.data=reshape(F.data,[F.width,F.height,F.depth]);
k=284+F.width*F.height*F.depth*fbyte;
[k,length(B),length(B)-k]
284byteまでがヘッダーです。Analyze formatみたいです。その後３次元画像データ部分になります。さらにその次があるのですが、よくわかりません。なくても大丈夫のようです。
ヘッダー部分、284byteあるんですが、ほとんど使用されてないので91〜284 byteはゼロになっています。
r,a,sは右、前、上の略です。
結果です。
F =
フィールドをもつ
struct:
version: 1
width: 256
height: 256
depth: 256
nframes: 1
type: 0
dof: 0
goodRASFlag: 1
spacingx: 1
spacingy: 1
spacingz: 1
xr: -1
xa: 0
xs: 0
yr: 0
ya: 0
ys: -1
zr: 0
za: 1
zs: 0
cr: 0
ca: 12.4000
cs: 19
ans =
1×194 の uint8 行ベクトル
1 列から 28 列
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
29 列から 56 列
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
57 列から 84 列
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
85 列から 112 列
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
113 列から 140 列
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
141 列から 168 列
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
169 列から 194 列
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
ans =
1×3 の int32 行ベクトル
16777500 16790218 12718
>>
読み込んだ画像を表示させてみます。
figure;
subplot(131);
imagesc(squeeze(F.data(round(F.width/2),:,:)));
daspect([F.spacingy,F.spacingz,F.spacingx]);
subplot(132);
imagesc(squeeze(F.data(:,round(F.height/2),:)));
daspect([F.spacingz,F.spacingx,F.spacingy]);
subplot(133);
imagesc(squeeze(F.data(:,:,round(F.depth/2))));
colormap(gray(256));
daspect([F.spacingx,F.spacingy,F.spacingz]);
結果です。
データ領域意向をナシにしたファイルを作成しました。
matlabでgzipで圧縮、gunzipで解凍となります。拡張子が.mgh.gzとなりますが、movefileでファイル名を.mgzに変更しています。
savename=[filename(1:end-4),'s.mgh'];
fid=fopen(savename,'wb');
fwrite(fid,B(1:(284+F.width*F.height*F.depth*fbyte)),'uint8');
fclose(fid);
gzip(savename);
movefile([savename,'.gz'],[savename(1:end-1),'z']);
delete(savename);
FreeSurfer6のfreeviewで見てみました。
FreeSurfer6のtkmeditで開いてみました。
データ領域部分以降がなくてもちゃんと認識されました。
surfのlh.pial
以下を実行します。
filename='c:\akira\surf\lh.pial';
fid=fopen(filename,'rb');
B=fread(fid,'*uint8');
fclose(fid);
k=1:3;
B(k)'
for kk=1:(length(B)-1)
if B(kk)==10
if
B(kk+1)==10
break;
end
end
end
kk=kk+1;
k=4:kk;
char(B(k))'
k=kk+(1:4);
n_vertices=swapbytes(typecast(B(k),'int32'))
k=kk+4+(1:4);
n_faces=swapbytes(typecast(B(k),'int32'))
k=kk+8+(1:(n_vertices*4*3));
vertices=swapbytes(typecast(B(k),'single'));
vertices=reshape(vertices,[3,n_vertices]);
k=kk+8+n_vertices*4*3+(1:(n_faces*4*3));
faces=swapbytes(typecast(B(k),'int32'));
faces=reshape(faces,[3,n_faces]);
[length(B),58+n_vertices*3*4+n_faces*3*4]
いろいろ試行錯誤して得た結果です。
big endianになっています。
最初はuint8 の3byte、255 255 254
次の4〜kk byteはテキスト文です。作者によってバイト数が変化します。
kkのところは10,10の改行コードです。
次はint32 頂点数、int32 面数です。
次はsingleで3×頂点数で頂点の座標、
次はint32で3×面数で三角メッシュの組み合わせです。0,1,2・・・の順番ですので、MATLABでは1を足して1,2,3・・・とします。
最後はテキスト文です。
ans =
1×3 の uint8 行ベクトル
255 255 254
ans =
'created by centos on Thu Jun
6 04:29:29 2013
'
n_vertices =
int32
160170
n_faces =
int32
320336
ans =
1×2 の int32 行ベクトル
5767430 5766130
左皮髄境界を表示します。
figure;
faces1=faces+1;
h=patch('vertices',vertices','faces',faces1');
set(h,'edgecolor','none','facecolor',[1,1,0]);
view([-120,20]);
light('position',[0,0,1]);
daspect([1,1,1]);
rotate3d on;
残りのテキスト文です。なくても問題ないようです。
k=(kk+8+n_vertices*3*4+n_faces*3*4+1):length(B);
char(B(k))'
結果です。
>>
ans =
' [1] valid
= 1 # volume info valid
filename = ../mri/filled-pretess255.mgz
volume = 256 256 256
voxelsize = 1.000000000000000e+00 1.000000000000000e+00 1.000000000000000e+00
xras = -1.000000000000000e+00
0.000000000000000e+00 0.000000000000000e+00
yras = 0.000000000000000e+00 0.000000000000000e+00 -1.000000000000000e+00
zras = 0.000000000000000e+00
1.000000000000000e+00 0.000000000000000e+00
cras = 0.000000000000000e+00
1.240000915527344e+01 1.900000000000000e+01
ƒmris_remove_intersection ../surf/lh.orig ../surf/lh.orig ProgramVersion: $Name: stable5 $ TimeStamp:
2013/06/05-17:37:26-GMT BuildTimeStamp: May 22 2011 08:23:31 CVS: $Id: mris_remove_intersection.c,v
1.6 2011/03/02 00:04:32 nicks Exp $
User: centos Machine: localhost.localdomain
Platform: Linux PlatformVersion: 2.6.32-71.el6.x86_64 CompilerName:
GCC CompilerVersion:
30400
•mris_make_surfaces -white NOWRITE -mgz -T1 brain.finalsurfs hashizume_akira lh ProgramVersion: $Name: stable5 $ TimeStamp:
2013/06/05-19:20:06-GMT BuildTimeStamp: May 22 2011 08:23:31 CVS: $Id: mris_make_surfaces.c,v
1.127 2011/03/02 00:04:33 nicks Exp $
User: centos Machine: localhost.localdomain
Platform: Linux PlatformVersion: 2.6.32-71.el6.x86_64 CompilerName:
GCC CompilerVersion:
30400 '
>>
なおGIfTIからmris_convertで変換したときは
' [1] valid
= 0 # volume info invalid
filename =
volume = 0 0 0
voxelsize = 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
xras = 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
yras = 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
zras = 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
cras = 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
'
となります。
テキスト文なしのファイルを作成してみます。
kk=1:(58+n_vertices*3*4+n_faces*3*4);
savename=[filename,'s'];
fid=fopen(savename,'wb');
fwrite(fid,B(kk),'uint8');
fclose(fid);
FreeSurfer6のfreeviewで読み込ませてみました。
FreeSurfer6のtksurferで読み込ませてみました。
テキスト部分がなくてもちゃんと開きました。
surfのlh.curv.pial
メッシュにはcurv曲率？があります。lh/rh.whiteとlh/rh.inflatedはlh/rh.curv、lh/rh.pialはlh/rh.curv.pialが対応しています。以下を実行します。
filename='c:\akira\surf\lh.curv.pial';
fid=fopen(filename,'rb');
B=fread(fid,'*uint8');
fclose(fid);
k=1:3;
B(k)'
k=3+(1:4);
n_vertices=swapbytes(typecast(B(k),'int32'))
k=7+(1:4);
n_faces=swapbytes(typecast(B(k),'int32'))
k=11+(1:4);
n_=swapbytes(typecast(B(k),'int32'))
k=15+(1:(n_vertices*4));
curv=swapbytes(typecast(B(k),'single'));
size(curv)
[length(B),k(end)]
いろいろ試行錯誤して得た結果です。
big endianになっています。
最初はuint8 の3byte、255 255 255
次はint32 頂点数、int32 面数、int32 何か？で、
最後がsingle 曲率×頂点数となっています。
ans =
1×3 の uint8 行ベクトル
255 255 255
n_vertices =
int32
160170
n_faces =
int32
320336
n_ =
int32
1
ans =
160170
1
ans =
1×2 の int32 行ベクトル
640695 640695
以下を実行します。
figure;
faces1=faces+1;
sulcus=curv;
h=trimesh(faces1',vertices(1,:),vertices(2,:),vertices(3,:),sulcus);
set(h,'facecolor','interp');
set(gca,'clim',[-1,1]*0.2);
colormap(jet(256));
daspect([1,1,1]);
view([-120,20]);
colorbar('veretical');
lh.inflatedとlh.curvを使って、mne_analyzeの画像を作成してみました。
filename='c:\akira\surf\lh.inflated';
fid=fopen(filename,'rb');
B=fread(fid,'*uint8');
fclose(fid);
for kk=1:(length(B)-1)
if B(kk)==10
if
B(kk+1)==10
break;
end
end
end
kk=kk+1;
k=kk+(1:4);
n_vertices=swapbytes(typecast(B(k),'int32'));
k=kk+4+(1:4);
n_faces=swapbytes(typecast(B(k),'int32'));
k=kk+8+(1:(n_vertices*4*3));
vertices=swapbytes(typecast(B(k),'single'));
vertices=reshape(vertices,[3,n_vertices]);
k=kk+8+n_vertices*4*3+(1:(n_faces*4*3));
faces=swapbytes(typecast(B(k),'int32'));
faces=reshape(faces,[3,n_faces]);
filename='c:\akira\surf\lh.curv';
fid=fopen(filename,'rb');
B=fread(fid,'*uint8');
fclose(fid);
k=3+(1:4);
n_vertices=swapbytes(typecast(B(k),'int32'));
k=7+(1:4);
n_faces=swapbytes(typecast(B(k),'int32'));
k=11+(1:4);
n_=swapbytes(typecast(B(k),'int32'));
k=15+(1:(n_vertices*4));
curv=swapbytes(typecast(B(k),'single'));
figure;
faces1=faces+1;
sulcus=curv;
h=trimesh(faces1',vertices(1,:),vertices(2,:),vertices(3,:),sulcus);
set(h,'facecolor','interp');
set(gca,'clim',[-1,1]*0.5);
map=flipud(gray(2)*0.5)+0.25;
colormap(map);
daspect([1,1,1]);
view([-120,20]);
colorbar('veretical');
このlh.curvやrh.curvがないとmne_analyzeでinflatedやwhiteを開くと
となります。lh.curv、rh.curvがあると
となります。pialの場合は、あんまり関係ないみたいです。lh.curvとrh.curv、lh.curv.pialとrh.curv.pialがないときです。
lh.curv.pialがrh.curv.pialがあるときです。違いはわかりません。
