提取ROI中的灰质体积

[转自52brain]

Good Ref 原文

鉴于有一些滴友询问如何提取特定ROI内的灰质体积，在这个帖子中我贴出来我使用的一个示例程序。供大家参考。这段程序其实基本和计算全脑灰质体积的程序一样，唯一不同的是使用了ROI，并且在计算灰质体积的同时读入ROI，并将ROI的每一层和灰质图的每一层点乘。

需要注意的是，ROI需要事先进行重采样（reslice），使它和你的灰质图维度一样。重采样的方法严超赣老师在别的地方讲过，这里就不详细说了。 好，废话不说，上干货：

clear;
subjIDs=load('I:\kongliliushuang\subjectID_list.txt'); %我把被试的编号都放在了一个txt文档里面（个人习惯而已，大家可以用别的方法。）
roi_dir = ('I:\kongliliushuang\ImagingData\ROI_img\'); % ROI所在的路径
cwd='I:\kongliliushuang\ImagingData/T1Img'; % 数据所在的路径
cd (cwd);
filename=dir ('smwc*.nii');     %我要处理的灰质文件都是smwc打头的，dir这个命令是列举文件夹下以smwc打头的文件

filestruct=[];
for i=1:size(filename,1)
   filestruct=[filestruct;{filename(i).name}]    %做一个smwc开头的文件名的列表，是cell形式的
end
    %%%%%%%%%%%%%%%%%%%%%Binary ROIs%%%%%%%%%%%%%%%%%%%%%
myROI={'insula.img','amygdala.img'};     % 重采样之后ROI图。可以根据需要放多个ROI

allVol=[];  %装最后体积结果的矩阵

for ii=1:size(myROI, 2)    %循环不同的ROI

ROI=[roi_dir,myROI{ii}];   % 到ROI 的路径
Vroi=spm_vol(ROI);  % 读进来ROI图


subj=[];       % 装每个ROI中所有被试结果的矩阵
allSubID=[];   %装被试编号的矩阵
for sub=1:size(subjIDs,1)    %循环所有被试
                tis=[];  %装每个被试全脑的体积的矩阵
                for tissue=1:1,  % tissue types, 1 是灰质
                        V=spm_vol(sprintf('c%d%d.nii',tissue,subjIDs(sub,1)))
                        tot=0;
                        for i=1:V(1).dim(3),   % V(1).dim(3)是slice的数量，下面这一段是计算每一个层（slice）的数据
                            img=spm_slice_vol(V(1),spm_matrix([0 0 i]),V(1).dim(1:2),0);    %读进来灰质图的特定的某一层
                            imgroi=spm_slice_vol(Vroi(1),spm_matrix([0 0 i]),Vroi(1).dim(1:2),0);   % 读进来ROI图的特定的某一层
                            new= img.*imgroi ;  % 两个层进行点乘，仅剩下ROI中非零的部分
                            tot=tot+nansum(new(:));   % 计算每一层中ROI中的数值，最后加起来
                        end;

                        voxvol=abs(det(V(1).mat))/100^3    %一个完完整整的体素的体积，升为单位
                        Vols=tot*voxvol    % 总体积                        
                        tis=[tis,Vols]
                end % tissue的结尾
           subj=[subj;sum(tis)];
           allSubID=[allSubID;subjIDs(sub,1)]
end   % subjects的结尾

allVol=[allVol,subj];
end   %ROI的结尾