Embedded Files
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のままで閲覧できます。
#######################################################
次に、ふたつのfeatures, A, Bに対応する後向推論画像のA AND NOT BとB AND NOT Aを取って画像出力すると言う課題を赤間研の学生が行ったものです(neurosynth_A_And_Not_B.py)。
ポイントは
or_ima = ima_first + ima_second
B_m_ima = or_ima + ima_second
M_b_ima = or_ima + ima_first
M_b_ima[M_b_ima>or_ima] = 0
B_m_ima[B_m_ima>or_ima] = 0
Pythonの場合、numpyのarrayであれば、MATLABと似た手が使えます。
Page updated
Google Sites
Report abuse