I. 2つの機能画像のAND, ORを取って新しい、画像を作成する方法
Neurosynthによる後向推論画像のAND領域, OR領域を計算する例でやってみよう。
Neurosynthから、emotion, rewardをfeaturesに、emotion_pFgA_z_FDR_0.01.nii.gzとreward_pFgA_z_FDR_0.01.nii.gzをダウンロード。
脳科学の世界で、この二つの語は類義語。当然情報が大きく重なっていることが予想される。その重なり具合を見て、双方に共通する(どちらか少なくとも一方に現れる)voxelsに関して、新たな画像を生成、閲覧、分析する。
PythonもMATLABも対話形式のインタープリタ言語であり、プログラム全体を一括してコンパイルすることなく、そのつどユーザーが、ソースコードを部分部分で逐次的に実行できる、使いやすい高水準のプログラミング言語です。これらの言語は特に配列の処理に適していますが、MRIの脳画像データも、デジタル化されている以上は多次元配列の形を取り、要素間(多くは対応するvoxels間)で、四則演算が可能です。以下のファイルでは、対応する要素間で値の積を取っており、出力した配列を脳画像データに戻しています。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I. MATLAB-SPMを使う方法。
Pythonと違ってMATLABはgzipで圧縮したgzファイルはそのままでは読み込めないので,
Ubuntu or Debian上の(Virtual Machineを含む)Terminalを使って解凍。
以下-dcオプションで元ファイルは残す。
gzip -dc emotion_pFgA_z_FDR_0.01.nii.gz >emotion_pFgA_z_FDR_0.01.nii
gzip -dc reward_pFgA_z_FDR_0.01.nii.gz >reward_pFgA_z_FDR_0.01.nii
これでniiファイルに解凍されている。
MATLABを起動し、その上で、spm fmriコマンドでSPM12を起動。
あらかじめaddpathコマンドなどで、spm12にパスを通しておくと良い。
今回はGUIでなくてMATLABのコマンドウィンドウ上で、SPMをコマンドラインで使う。
以下
SPM(Statistical Parametric Mapping)関係
1)SPMの一口知識(少々古いです。SPMをGUIからでなく、batch処理用のscriptを書きたい人向け)
も参照する。
https://sites.google.com/site/akamatitechlab/research/tips/spm1
[この質問を参考にしてください。
(質問)SPMでコマンドラインを用い、あるいは自作のスクリプトを書いて、fMRI画像中のvoxelsの賦活値を抜き出したり、条件に合ったvoxelsのindex情報を取得したり、賦活値を変更したりすることはできますか？
(回答)以下のようにすれば簡単にできます。
V=spm_vol('filename.img');%Analyze(あるいはNifti)ファイルを読み込む
ima=spm_read_vols(V);%各voelsの値を抜き出す。返ってくるのは3次元の配列
ima(ima>threshold)=1;
ima(ima<=threshold)=0;
%thresholdより大きい値は1、以下の値は0になります。
V.fname='newfilename.img';
spm_write_vol(V,ima);%新しい画像に保存します。]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
datalist={'emotion','reward'}%セル配列
%https://jp.mathworks.com/help/matlab/cell-arrays.html
%セル配列の概念を理解する。
first_fname=strcat(datalist{1},'_pFgA_z_FDR_0.01.nii')%文字列の水平配列;niiファイル名の文字列
%emotion_pFgA_z_FDR_0.01.nii
second_fname=strcat(datalist{2},'_pFgA_z_FDR_0.01.nii')%文字列の水平配列
reward_pFgA_z_FDR_0.01.nii
and_fname= strcat('and_',datalist{1},'_',datalist{2},'.nii')%文字列の水平配列;出力ファイル名を準備
%and_emotion_reward.nii
or_fname= strcat('or_',datalist{1},'_',datalist{2},'.nii')%文字列の水平配列
%or__emotion_reward.nii
first_V=spm_vol(first_fname);%Niftiファイルを読み込む
second_V=spm_vol(second_fname);%Niftiファイルを読み込む
%;(セミコロン)を付けずにfirst_Vが何を返すか確認
ima_first=spm_read_vols(first_V);%各voelsの値を抜き出す。返ってくるのは3次元の配列
ima_second=spm_read_vols(second_V);%各voelsの値を抜き出す。返ってくるのは3次元の配列
and_ima=ima_first.*ima_second;
%2つの配列の乗算.https://jp.mathworks.com/help/matlab/ref/times.html
%要素ごとの積なので、どちらかひとつが0だと結果は0でANDが取れる。ふたつのfeaturesでともに有意なvoxelsのz-scoreの積が計算できる。
and_ima(and_ima>0)
%で0より大きい値が検索できる。配列で条件に基づいて値を抜き出すやり方
%size(and_ima(and_ima>0),1)が946を返すので946のvoxelsが共通。
or_ima=ima_first + ima_second;
%こちらは要素ごとの和
%size(or_ima(or_ima>0),1)が 13476を返すので、13476個のvoxelsがどちらか少なくとも一方で0より大きい値を持っていることが分かる。
first_V.fname=and_fname;
%これにより、spm_vol関数が返す構造体から参照されるファイル名が、ANDの計算結果を代表するものに更新される。
%言いかえると、voxelの計算結果とファイル名以外のメタデータ情報は、first_Vそのままが継承できる。
spm_write_vol(first_V, and_ima);
%ANDの結果and_imaは新しい画像'and_emotion_reward.nii'に保存されます。
first_V.fname=or_fname;
%これにより、spm_vol関数が返す構造体から参照されるファイル名が、ORの計算結果を代表するものに更新される。
%言いかえると、voxelの計算結果とファイル名以外のメタデータ情報は、first_Vそのままが継承できる。
spm_write_vol(first_V, or_ima);
%ORの結果or_imaは新しい画像'or_emotion_reward.nii'に保存されます。
%Displayボタンから新しい画像を見てみよう。
%XjViewから新しい画像を見てみよう。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
###########################################################
II. Pythonとその様々なパッケージを用いる方法。
LinuxのTerminalからipythonを起動し、
以下の####と####の間の行を全部、クリップボードにコピーして
%paste
というマジックを利用します。
####################################################
import numpy as np
from nipy import load_image, save_image #Needs nipy
from nipy.core.api import Image
import os
# -*- coding: utf-8 -*-
datalist=['emotion','reward']
and_fname='and_' + datalist[0] + '_' + datalist[1] + '.nii.gz'
or_fname='or_' + datalist[0] + '_' + datalist[1] + '.nii.gz'
#Pythonは0から始まる。 + で文字列の連結ができる。
first_V = load_image(datalist[0] + '_pFgA_z_FDR_0.01.nii.gz')
ima_first = first_V.get_data()
second_V = load_image(datalist[1] + '_pFgA_z_FDR_0.01.nii.gz')
ima_second = second_V.get_data()
#この辺はspm_vol(), spm_read_volsとパラレルと見なして良いでしょう。
and_ima = ima_first * ima_second
or_ima = ima_first + ima_second
save_image(Image(and_ima, first_V.coordmap), and_fname)
save_image(Image(or_ima, first_V.coordmap), or_fname)
#######################################################
これをipython上でniiファイルに解凍したい場合は
os.system('gzip -cd '+ and_fname + ' > ' + '' + fname[:-3])
ただし、FSLViewなどでは、nii.gzのままで閲覧できます。