NIRS_SPM 批处理成像的一个小改进

This is a guest post by Dr. Ning Liu from Stanford.

前些天 在用 NIRS_SPM批处理 做fNIRS脑图时，有几组数据总是报错（大约占总数的三分之一左右），总是说无法生成图像。

??? Error using ==> image
Error using ==> image
Image CData can not be complex

Error in ==> imagesc at 19
hh = image(varargin{1},'CDataMapping','scaled');

Error in ==> activation_map_batch at 365
imagesc(stat_brain);

开始以为是前面几步处理时出了什么问题，于是认认真真从头又做一遍，问题照旧。焦头烂额几天，发现做插值(interpolation)之前一切都是美好的，问题就是在插值以后生成了 复数的interp_var！没办法，只好调出NIRS_SPM 里的批处理函数(activative_map_batch.m)，一步一步跟进追踪。Activative_map_batch.m里面 是这么计算interp_var的：

nch = length(chs{kk}); % # of channels
mtx_var = zeros(nch);
for aa = 1:nch
for bb = 1:nch
mtx_var(aa,bb) = var(chs{kk}(aa), chs{kk}(bb));
end
end
[V_X D_X] = eig(mtx_var);
tmp = D_X.^(1/2) * V_X' * B{kk};
interp_var = [interp_var sum(tmp.^2,1)];

问题就出在eig函数那里。Eig函数是用来求本征值和本征向量的，但是在运行某些数据时，会生成一些很小很小的本征值，通常远小于别的本征值十个数量级以上，非常接近于零。我的情况下生成的本征值就是1e-20 量级。 这些奇怪的值（特别是负数的极小值）导致了下面一步开平方出现了复数(complex number),从而报错提示不能使用图形生成函数imagesc。这些值的产生通常是因为做数字近似的时候产生的误差(artifact of numerical approximation)。找到原因，解决办法就很简单了，把这些极小值用零替换就好了。改进后的程序是这样的：

nch = length(chs{kk}); % # of channels
mtx_var = zeros(nch);
for aa = 1:nch
for bb = 1:nch
mtx_var(aa,bb) = var(chs{kk}(aa), chs{kk}(bb));
end
end
[V_X D_X] = eig(mtx_var);

% round very small engenvalues to zero <--by NL
% very small engenvalues result from eig function
% usually is an artifact of numerical approximation. simply
% round them to zero here.
[D_X_err_I1, D_X_err_J1] = find(D_X <0);
disp(D_X(D_X_err_I1, D_X_err_J1));
D_X(D_X_err_I1, D_X_err_J1) = 0;
[D_X_err_I2, D_X_err_J2] = find(D_X>0 & D_X<1e-15);
disp(D_X(D_X_err_I2, D_X_err_J2));
D_X(D_X_err_I2, D_X_err_J2) = 0;
% #

tmp = D_X.^(1/2) * V_X' * B{kk};
interp_var = [interp_var sum(tmp.^2,1)];

本文作者：刘宁博士

如果需要改进后的批程序，请发电子邮件给本文作者刘宁博士: [email protected]

---

---

写作助手，把中式英语变成专业英文

---

Want to receive new post notification? 有新文章通知我