查看: 43858|回复: 146
打印 上一主题 下一主题

[spm操作] VBM分析中,如何计算全脑的灰质、白质和脑脊液的体积 ? [复制链接]

Rank: 7Rank: 7Rank: 7

水晶
6004
心级
2643
精华
4
主题
36
帖子
2242
跳转到指定楼层
楼主
发表于 2014-12-3 04:57:45 |只看该作者 |倒序浏览
本帖最后由 空里流霜 于 2019-7-31 09:55 编辑

本帖作为《用Matlab和SPM批量处理被试的经验总结》的一部分
目录贴请见http://home.52brain.com/forum.ph ... =1&extra=#pid158525

在做VBM分析的时候,经常会把全脑体积(total brain volume=灰质+白质)或是颅内体积(total intracranial volume=灰质+白质+脑脊液)作为无关变量放在第二层的GLM中。
因此计算灰质白质脑脊液的体积就成了一个有用的步骤。
我参照John Ashburner在SPM的mailing list中做的一小段程序,稍作修改(几乎没改)放在此处供大家参考。至于在SPM中如何做VBM,大家可以参考此教程:www.fil.ion.ucl.ac.uk/~john/misc/VBMclass10.pdf

如果你用VBM8这一工具包的话,里面会同时生成一个文本文档,里面就有被试的这三种组织的体积。但VBM8估算出的体积和Ashburner的代码估计出的有一些出入,具体可能是计算方法的不同。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%此程序来自于Ashburner在SPM mailing list中的回复,我做了些许修改以适用于我的数据以及批处理的要求
clear;
cd('D:\VBM\Dartel\results')  %数据所在之处
subj=[];
for sub=1:20??%多少个被试
tis=[];
for tissue=1:3,    %三种组织,灰质,白质,脑脊液
? ?? ???V=spm_vol(sprintf('c%dsubject%02d_T1.nii',tissue,sub))  %可以是nii格式,也可以是img格式
? ?? ???Vols=zeros(numel(V),1);
? ?? ???for j=1:numel(V),
? ?? ???tot=0;
? ?? ???for i=1:V(1).dim(3),
? ?? ?? ?? ?img=spm_slice_vol(V(1),spm_matrix([0 0 i]),V(1).dim(1:2),0);
? ?? ?? ?? ?img=img(isfinite(img)); %exclude non-finite values
? ?? ?? ?? ?tot=tot+sum(img(:));
? ?? ???end;
? ?? ???voxvol=abs(det(V(1).mat))/100^3 %一个voxel的体积,以升为单位,必须加上这一句
? ?? ???tot=tot %integral of voxel intensities
? ?? ???Vols(j)=tot*voxvol
? ?? ???end
? ?? ???tis=[tis,Vols]
end
subj=[subj;tis]
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
在这里我说一下上面程序的思路。你会发现有脑组织的体素的intensity是1或者接近1.如下图就是第一百层的轴状位的大脑(右侧),虽然它是个图,但本质上是个矩阵(左侧)。你会看到很多地方是0,那是非脑区的地方,也有很多地方是接近1的。矩阵里面向左凹进去那一块其实就是大脑额叶前面凹进去的那部分。
你用这一句也能看到第一百层的矩阵 img=spm_slice_vol(V(1),spm_matrix([0 0 100]),V(1).dim(1:2),0);

Ashburner可能假设如果一个体素是1,那么脑组织就占这个体素的百分之百,如果这个体素的intensity是0.5,那么脑组织就占百分之五十。所以把所有的体素加起来就是大脑中的以体素为单位的体积。
但体素的大小在每个研究中可能都是不太一样的。尤其是使用c1开头的刚分割完的比较原始的图像时,比如我这个图就是体素大小0.479×0.479×1。所以还要再通过体素的大小换算成具体的毫升或是升等单位。
而这一句voxvol=abs(det(V(1).mat))/100^3就是做这个事的。V(1)是我们载入的图像本身,在这里是个结构体(structure),而V(1).mat是一个仿射转换矩阵(affine transformation matrix),里面记录了体素的大小的信息。det函数做的就是将体素的体积信息提取出来。abs就是看绝对值的,因为det提取出的值可能是负的。







有的同学会问“我是用标准化之前的图像还是标准化之后的图像来计算这个体积?”或者“我是用mudulation之前的还是之后的?”或者“我是用smooth之前的还是之后的”。我用自己的数据算了一下,其实我感觉是这样的。用smwc1开头的文件(经过了标准化,modulation和smooth)和用c1开头的文件(刚分割完的原始空间的灰质)来计算的话,灰质的绝对体积值是不一样的。但是被试之间的相对体积是非常接近的,四舍五入到小数点后三位数是一样的(由于smooth肯定影响了具体的数值,smooth之前应该更加接近)。比如,有五个被试用c1开头的文件算出来的数值依次是
0.6083
0.8447
0.6546
0.7670
0.8260

用smwc1开头的文件算出来的数值依次是
??0.6053
??0.8408
??0.6515
??0.7621
??0.8217

可以看出来还是很接近的。
另外对于同样的5名被试,我还看了VBM8处理后生成的txt文件里面的灰质体积,虽然感觉有点差别,但体积还是基本一致的,
594.832
823.331
633.403
746.745
754.158

我还没有尝试标准化modulation但未smooth的情况。所以我想标准化和modulation并没有保存其绝对体积值,而是保留了相对体积值。

如果你要把体积作为无关变量放在GLM模型中的话用标准化和调制之后的是没有问题的,但要计算绝对体积的话,那个c1开头的图像,就是那个在原始空间的刚刚分割完的图像算出来的体积是最靠谱的。

【代码中sprintf函数的用法】
这里附带说明一下sprintf函数的用法,举个例子说sprintf('c%dsubject%02d_T1.nii',tissue,sub)中,红色的地方是%d,是说%d的地方由后面红色的部分(tissue)的变量所代替,
蓝色的部分是%02d,是说%02d的地方由后面蓝色部分(sub)的变量所代替,其中的02是说把sub变量转换成两位数,比如sub是5,转换之后就是05,如果sub是13,转换之后还是13。所以你能看出来,上面的%是表示说这里要填上后面某个对应的变量,d表示是整数,02是表示要保持数字是几位的。
这样,如果sub是4,tissue是3,那么这个sprintf('c%dsubject%02d_T1.nii',tissue,sub)输出之后就是c3subject04_T1.nii,而这个就是我的数据的命名方式。
如果你的命名是c1s-0003-00001-000001-01.nii,需要把组织对应的1换成%d,被试对应的部分,比如0003换成%04d。
每一部分每一个符号都是有用的,需要你仔细研究一下。也可以在matlab的命令窗口自己先尝试一下,这样方便debug。比如你在命令窗口输入 sprintf('%02d',5),看会输出什么。如果换成 sprintf('%03d',15),又会输出什么。这样你就能有一点感觉这个函数是如何使用的。







附件: 你需要登录才可以下载或查看附件。没有帐号?注册
已有 3 人评分水晶 收起 理由
xhhyjm11 + 1 赞一个!
neal仰望 + 1 很给力!
edengao007 + 1 很给力!

总评分:?水晶 + 3? ?查看全部评分

转发到微博
空里流霜不觉飞,汀上白沙看不见。

Rank: 2

水晶
157
心级
35
精华
0
主题
3
帖子
27
沙发
发表于 2014-12-3 16:16:42 |只看该作者
大神,请问计算过后三种物质的体积是以什么为单位的?还有,对计算体积所用到的数据是否有要求,比如必须在做smooth之前才能进行体积计算等等?我找到两个脚本,但是计算的结果不同,请问这是怎么回事?谢谢!

Rank: 7Rank: 7Rank: 7

水晶
6004
心级
2643
精华
4
主题
36
帖子
2242
板凳
发表于 2014-12-4 09:02:28 |只看该作者
我不爱脑科学 发表于 2014-12-3 03:16
大神,请问计算过后三种物质的体积是以什么为单位的?还有,对计算体积所用到的数据是否有要求,比如必须在 ...

我上面这一段是以litre为单位的。
我也遇到过类似的问题,我是在标准化之前和之后计算的,好像是不太一样。我个人觉得在标准化之前计算的应该最接近真正的体积。
空里流霜不觉飞,汀上白沙看不见。

Rank: 3Rank: 3

水晶
110
心级
82
精华
0
主题
18
帖子
46
地板
发表于 2014-12-4 09:33:31 |只看该作者
空里流霜 发表于 2014-12-4 09:02
我上面这一段是以litre为单位的。
我也遇到过类似的问题,我是在标准化之前和之后计算的,好像是不太一样 ...

我用get-total获得的是ml为单位的,请问师兄ml和体素之间怎么换算啊?还有用VBM做的话,生成的那文本文档里得到的灰质白质脑脊液的值得单位也是ml吗?

Rank: 7Rank: 7Rank: 7

水晶
6004
心级
2643
精华
4
主题
36
帖子
2242
5
发表于 2014-12-6 01:00:18 |只看该作者
本帖最后由 空里流霜 于 2014-12-11 22:25 编辑
D2012LY 发表于 2014-12-3 20:33
我用get-total获得的是ml为单位的,请问师兄ml和体素之间怎么换算啊?还有用VBM做的话,生成的那文本文档 ...

我原来说得不对,现在改正过来。在上面的代码中,abs(det(V(1).mat))其实算出来的是体素的以立方毫米为单位的体积,为了转换成升,就要除以100的三次方。若想以毫升为单位就要把100的三次方改成10的三次方。
空里流霜不觉飞,汀上白沙看不见。

Rank: 3Rank: 3

水晶
110
心级
82
精华
0
主题
18
帖子
46
6
发表于 2014-12-8 08:40:49 |只看该作者
空里流霜 发表于 2014-12-6 01:00
这个要看你的体素大小。我上面的程序是认为你的体素是1*1*1立方毫米的。所以如果要以毫升为单位,每个值 ...

这是要怎么知道体素的值得啊?在xjview里的体素是什么单位的师兄知道吗?

Rank: 7Rank: 7Rank: 7

水晶
6004
心级
2643
精华
4
主题
36
帖子
2242
7
发表于 2014-12-9 11:00:55 |只看该作者
D2012LY 发表于 2014-12-7 19:40
这是要怎么知道体素的值得啊?在xjview里的体素是什么单位的师兄知道吗?
...

用SPM界面的Display按钮就可以看到。在右下方会有vox size的字样。后面的值就是体素大小。
空里流霜不觉飞,汀上白沙看不见。

Rank: 3Rank: 3

水晶
266
心级
111
精华
0
主题
2
帖子
93
8
发表于 2014-12-9 16:09:04 |只看该作者
matlab 高手,佩服,俺们matlab就是不懂,好好学习。非常感谢。

Rank: 7Rank: 7Rank: 7

水晶
6004
心级
2643
精华
4
主题
36
帖子
2242
9
发表于 2014-12-11 11:33:00 |只看该作者
xmzshljz321 发表于 2014-12-9 03:09
matlab 高手,佩服,俺们matlab就是不懂,好好学习。非常感谢。

客气!多交流!
空里流霜不觉飞,汀上白沙看不见。

Rank: 3Rank: 3

水晶
110
心级
82
精华
0
主题
18
帖子
46
10
发表于 2014-12-12 11:15:12 |只看该作者
空里流霜 发表于 2014-12-9 11:00
用SPM界面的Display按钮就可以看到。在右下方会有vox size的字样。后面的值就是体素大小。
...

这下彻底对上了,谢谢师兄
您需要登录后才可以回帖 登录 | 注册

bottom

Powered by Discuz! X2

? 2001-2011 Template By Yeei. Comsenz Inc.

回顶部