VBMで得られる画像を用いて、灰白質や白質画像の容積を求めることができます。
事前の準備として、以下のことが必要です。
- get_totals.m
- WFU PickAtlas
Ged Ridgwayが公開している容積算出のためのスクリプトです。http://www0.cs.ucl.ac.uk/staff/g.ridgway/vbm/get_totals.mここから入手できます（右クリックしてリンク先を保存としてください）。これをSPMのディレクトリに保存します。
これは、局所容積を求めたい時に便利です。マスク画像を簡単に作成することができます。http://fmri.wfubmc.edu/software/PickAtlasここから入手できます。今回はこのセットアップ方法は説明しませんのであしからず。
画像としては、以下のものが必要となります。
- SPM8のNew Segmentを使っている時
- SPM12のNew Segmentを使っている時
- VBM8を使っている時
modulationをかけ、smoothを[0 0 0]としたもの。出力ファイルは、灰白質ならば、smwc1xxx.niiのような形をとるはずです。New Segmentでは、smoothingをかけないという選択肢はないので、[0 0 0]とします。今までいろいろ試してみましたが、smoothingをかけると信号値が若干変化しますので、smoothing前のデータを用いるのがよいと思います。
modulationをかけ、smoothを[0 0 0]としたもの。出力ファイルは、灰白質ならば、mwc1xxx.niiのような形をとるはずです。SPM8との違いは、smoothingを[0 0 0]とすると、sがつかないようになったようです。（猪狩先生、情報をどうもありがとうございます）
Modulationを”affine + non-linear (SPM8 default)”とします。non-linear onlyですと、線形変換の結果は考慮されませんので、頭蓋容積を補正済みのデータになってしまい、生値を得ることができないので注意が必要です。出力ファイルは、灰白質ならば、mwrp1xxx.niiとなります。
それでは、はじめましょう。
- 基本編
- 応用編
- 作業ディレクトリの準備
- マスク画像の準備
- マスク画像のボクセルサイズの変更
- get_totalsで左半球の灰白質容積算出
まずは、全灰白質容積を算出しましょう。
Matlabのコマンドウィンドウから、
get_totals
とタイプします。そうすると、画像を選択するダイアログが出ますので、灰白質容積を求めたい画像を選択します。（複数選択できます。）
すると、Matlabのコマンドウィンドウに、
ans = 739.8652
というように出力されます。これが求めたい容積です。複数画像を選択した場合、ファイルの選択順に出力されます。
これは簡単ですね。get_totalsとコマンドを打てばいいのですから。
それでは、次にいきましょう。
get_totalsでは、マスク画像を指定することにより、局所の容積を算出することができます。
ここでは、左半球の灰白質容積を求めたいとします。
Matlabのコマンドウィンドウで、
help get_totals
とタイプしてみてください。そうすると、次のようなメッセージが出てきます。
get_totals – Returns image totals (sum over all voxels), in ml
t = get_totals
[t files] = get_totals(files, thr, msk)
GUI file-selection is used if files not specified as argument (or empty).
ここにあるget_totals(files, thr, msk)がキモでして、get_totalsの後に閾値やマスク画像を指定できます。これを手で入力するのは面倒くさいのですが、[]と指定して、「空ですよ」と伝えてあげるとファイルの選択画像が出現します。なので、これを使いましょう。
それでは、具体的に進めていきます。
SPMで仕事をするときには、作業ディレクトリ(Working directory)を決めておいた方が便利です。Matlabのコマンドウィンドウなり、SPMのUtil->CDで、そのディレクトリにうつっておいて後の仕事をします。
WFU PickAtlasでマスク画像を作成します。SPM8のToolboxからWFU PickAtlasを起動し、Human Atlasを選択すると、下図のような画面が起動します。
今は、左半球のマスクを作りたいと考えています。そこで、PickAtlasの左側にあるメニューから、TD Hemispheresを選択し、Left Brainstem, Left Cerebellum, Left Cereberumをそれぞれダブルクリックして、右のパネル(Working Region1)に移します。「大脳の情報だけわかればいいから」という人はLeft Cereberumだけでもいいかもしれません。今は、複数の領域からマスクを作れることを示すためにあえて左を全部選びました。
そうしましたら、画面中央下あたりにある、”SAVE MASK”をクリックします。そうすると保存するための画面が出ますので、作業ディレクトリの中に適当に名前をつけて保存してください。ここでは、lt_hemisphere.niiという名前にします。そうすると、マスク画像が作成されます。
試しに、SPMのDisplayからlt_hemisphere.niiを表示してみましょう。下図のようになるはずです。ここで、Voxel Sizeが2x2x2となっていることに注意してください。
先ほど、Displayの結果で見たように、PickAtlasで作成されるマスク画像はボクセルサイズが2x2x2となっています。しかし、DARTELの出力結果のボクセルサイズは標準では、1.5×1.5×1.5となっています。これでは画像の次元が異なってしまうため、画像の掛け算を行うことができません。したがって、マスク画像のボクセルサイズを変更することが必要となります。このために、SPMの”Coregister (Reslice)”を用います。
SPMの左上のウィンドウから、”Coregister (Reslice)”を選択します。
そうすると、下図のようなウィンドウが出現しますので、
Image Defining Spaceに、灰白質画像を指定します。
Images to Resliceに、マスク画像を指定します。
InterpolationはNearest neighbourにします。そうすることで、マスクがバイナリ画像（二値画像）のままとなります。
その他の設定はデフォルトのままで結構です。
そうすると、計算が始まり、新しくファイルの先頭にrが付加されたファイルが生成されます。今の場合ですと、rlt_hemisphere.niiという画像ができます。
これをSPMのDisplayで表示してみます。そうすうと、Vox Sizeが-1.5 x 1.5 x 1.5と変わっていることがわかるかと思います。
ここまで来ればあと一息です。
Matlabのコマンドウィンドウで、次のようにタイプします。
get_totals([],[],[])
ここで、[]はMatlabではそこが空ということを意味します。今は、画像ファイルも閾値もマスク画像も空ですよという意味です。そうすると、最初に下図のように画像ファイルを選択するダイアログが表示されます。
そうすると、次にマスク画像を選択する画面が出てきます。ここでは、先程準備したrlt_hemisphere.niiを選択します。
なお、閾値は空ですと、0に設定されます。もし、信号値0.1以上を示す領域での容積を求めたいならば、ここで[]の代わりに、0.1と入力すれば大丈夫です。
そうすると、Matlab上に
ans = 350.6627
のような表示が出てきます。先ほどの全灰白質容積と比較すると約半分になっていることがわかりますね。
このようにしてVBMで得られた画像から容積を求めることができます。
いつも本サイトや先生の著書｢すぐできるVBM｣を拝読し勉強しております。
名古屋市立大学 脳神経内科の谷口と申します。
SPM12を用いてMRI画像のVBMを行っており、局所容積の計算や関心領域の評価に進もうとしていますが、WFU_PickatlasをダウンロードしSPMのtoolbox内でzipファイルを解凍してMatlabのパス設定で指定しても、WFUがtoolbox内に格納されず使用できない状態となってしまいます。
また根本先生が本ページにリンク先に掲載していただいているWFUのサイトも、インターネットセキュリティの関係なのか閲覧できない状態です。。
根本先生にこのような初歩的な段階で相談するのは大変僭越と存じますが、お手すきの際にご意見頂戴できると幸いです。
何卒よろしくお願い申し上げます。
谷口先生
ご質問ありがとうございます。
おっしゃるようにWFUのサイトはなくなっていますね。
今、NITRCからしかダウンロードできないようです。
https://www.nitrc.org/projects/wfu_pickatlas/
こちらから、
WFU_PickAtlas_3.zip
をダウンロードいただき、展開すると、
wfu_pickatlas
wfu_results
wfu_tbx_common
の3つのフォルダが入っていることが確認できるかと思います。
この3つのフォルダを、
SPM12 の toolbox 内に置いていただけますか。
Matlabのパス設定は追加しなくて大丈夫です。
それで、SPMを起動して、toolbox を選ぶと wfupickatlas が 出てくるのではないかと思います。
ご確認いただけますか。
根本清貴
根本先生
早速のご回答ありがとうございます
ご指示通り作業を行いましたが、spmのtoolbox内にwfupickatlasは読み込みされませんでした。
同じ方法でAAL3はtoolbox内に利用できました。考えられるトラブルは何かありますでしょうか？
たびたびの初歩的な質問となり恐縮です。
谷口葉子
谷口先生
– そうしましたら、まず、spm12/toolbox の直下に、wfu_pickatlas があり、その中に、wfu_pickatlas.m があるかどうかを確認してください。
– そのうえで、SPMを終了し、Matlabも終了し、改めて、Matlabを起動し、SPMを起動して確認してください。
それでいかがでしょうか？
もし、それでもだめならば、Matlabのパス設定で、spm12/toolbox/wfu_pickatlas を、spm12のパスの下になるように設定してみてください。
そのうえで、Matlabのコマンドウィンドウから
とタイプしたら起動するかと思いますがいかがでしょうか？
根本先生
教えて頂いた方法でwfupickatlasをtoolbox内に展開できました。
ご指導くださり、誠にありがとうございます。
谷口葉子
谷口先生
無事解決したようでよかったです。
いつも本サイトや根本先生の著書｢すぐできるVBM｣で勉強をさせて頂いております。
豊川市民病院の打田と申します。定量的磁化率画像（QSM）を用いてパーキンソン病の画像解析などをしています。
get_totalsはモジュレーション前の灰白質画像や白質画像のvolumeも測定できるのでしょうか。測定できる場合、QSMのvolumeもQSM.nii（QSMのNIfTI画像）をinputするだけで測定できるということでよろしいでしょうか。
また、例えばQSM.niiで画素値が0.1より高い箇所のvolumeを求める場合、get_totals([],0.1,[]) で解析すればよろしいでしょうか。この場合、マスク画像は何を選択すればよろしいでしょうか。
質問が重なり申し訳ありませんが、どうぞよろしくお願いいたします。
打田先生
すみません、コメントに気づくのが遅くなってしまいました。
1. get_totalsはモジュレーション前の灰白質画像や白質画像のボクセル値の合計も測定できます。
2. QSMに関して、私は直接触ったことがないのですが、ボクセル値=定量値であるならば、使えるはずです。
3. 画素値が0.1より高いかしょのvolumeを求めたい場合は、get_totals([],0.1,[])で大丈夫です。マスク画像については、ご自身で決める必要があります。通常は脳をカバーするマスク画像を使うことが多いかなと思います。
以上、よろしくお願いします。
根本清貴
根本先生
ご返信ありがとうございます．２に関して，ボクセル値＝定量値の場合，get_totalsから求まる値は，volumeではなくて”定量値の合計”という理解でよろしいでしょうか．この場合，”volume”や”定量値の平均”を求める方法はございますでしょうか．
打田佑人
打田先生
はい、そのとおりです。
ボクセル値が定量値の場合、Volumeは、ボクセル数×1ボクセルの体積 ということになります。
たとえば、ある領域が 3000 ボクセルあって、 1ボクセルが 2mm x 2mm x 2mm で構成されているならば、
体積は
3000 x 2 x 2 x 2 = 24000mm^3 = 24mL となるわけです。
それでは、ボクセル数をどう求めるかということになります。
これはマスク画像を使って求めることになります。
https://www.nemotos.net/?p=113
こちらにあるスクリプトをお使いください。
Matlabに1行ずつ入力すると動きが学べると思います。
file = spm_select;
ここでマスク画像を選んでください。
vol = spm_vol(file);
マスク画像のヘッダー情報を、Matlabで vol という変数に読み込みます。
img = spm_read_vols(vol);
マスク画像のヘッダー情報を用いて、Matlabで画像を直接 img という変数に読み込みます。
3次元の行列になります。
sum(img(:)>0)
img変数で、各行列で、値が0より大きいセルの数=ボクセルの数を数えます。
この結果がボクセル数になります。
根本先生
ありがとうございます．大変よくわかりました．追加ですみません，先のメールの1についてもご質問させてください．’1. get_totalsはモジュレーション前の灰白質画像や白質画像のvolumeも測定できる’ ということですが，見出しのテーマの事前の画像準備の段階では，modulation後の画像を用意する必要があるのは，なぜなのでしょうか．
打田佑人
打田先生
すみません、表現が正確ではありませんでした。
modulation前の画像のボクセル値は 「確率」を示します。
ボクセル値 0.8 は、80%の確率で灰白質があるという意味です。
一方、modulation後の画像のボクセル値は「容積」です。
ボクセル値 0.8 は、そのボクセルの80%を示す値が体積という意味になります。
1mm x 1mm x 1mm のボクセルであれば、0.8mm^3 ということになります。
なので、モジュレーション前の灰白質画像や白質画像のvolumeも測定できるというのは正確に言えば、
モジュレーション前の灰白質画像や白質画像のボクセル値の合計も計算できる
という意味になります。（先のコメントも修正しました）
根本清貴
根本先生
大変よくわかりました、質問が重なり申し訳ありません。https://www.nemotos.net/?p=798では、モジュレーション前のc1*.niiで灰白質容積を求める、と記載してございますが、こちらもモジュレーション後の方が望ましいのでしょうか。
打田佑人
打田先生
標準化する前と後で話が変わります。
標準化する前：c1*.nii は確率マップではあるのですが、標準化する前はほとんど容積を反映すると考えられているので、これを使ってもらって大丈夫です。
標準化した後：標準化した後は、確率＝容積では全くなくなってしまいます。標準化のプロセスで大きさが変わってしまうので。そのためにモジュレーションが必要になります。
なので、ご指摘いただいた記事に関しては、標準化前の情報を使うので、c1を使って大丈夫ということになります。
根本清貴
根本先生
ありがとうございます．いつも大変わかりやすい説明で感謝しております．今後ともどうぞよろしくお願い致します．
打田佑人
根本先生
たびたびすみません，先の質問の中で，また疑問点が出てきました．
“modulation後の画像のボクセル値は「容積」です。ボクセル値 0.8 は、そのボクセルの80%を示す値が体積という意味になります。1mm x 1mm x 1mm のボクセルであれば、0.8mm^3 ということになります。”とご教示頂きました．DARTELのボクセルサイズは標準で1.5×1.5×1.5ですから，DARTELで標準化およびmodulation後の画像のボクセル値＝0.8の場合，1.5×1.5×1.5×0.8mm^3が実際のボクセルの体積ということになりますでしょうか．そのように考えると，DARTELで標準化およびmodulation後の画像の容積を知りたい場合，get_totalsで出力されるボクセル値の合計に，ボクセルサイズ（1.5×1.5×1.5）を乗じる必要があるように思うのですが，いかがでしょうか．
打田佑人
打田先生
実際のボクセルの体積はおっしゃるとおりです。
そして、get_totalsは、1ボクセルの容積を勘案してくれています。
なので、構造画像においては、心配はいらないです。
（先生の例で言えば、1.5×1.5×1.5×0.8の値を出力してくれます）
定量画像などの場合は、それぞれ慎重に検討する必要があると思います。
根本清貴
根本先生
ありがとうございます．最後にすみません，先生は，「体積」と「容積」をどのように使い分けておられますでしょうか．
打田佑人
打田先生
すみません、使い分けていません…。
その場で頭に浮かんできた方を使っているだけなので気にされないでください。
根本清貴
根本先生
細かい質問になり申し訳ありません．根本先生の説明はいつも大変わかりやすくて助けられています．今後ともどうぞよろしくお願いいたします．
打田佑人
いつも先生のサイトを読ませて頂き、研究を行なっております。
東京医科歯科大 耳鼻科 山本と申します。
VBMで得られた画像より、get_totallsを使用したりしてROIの体積などは解析しております。
VBMでは密度も算出できるとなっておりますが、具体的にどの数値になるのかご教授頂けますでしょうか？
得られた体積を使用して、算出するのでしょうか？
基本的な質問で申し訳ないことでございます。
山本先生
すみません、先生のコメントがなぜかスパム扱いされており、今しがた気づきました。
VBMでの密度(density)は、純粋に0-1の間となります。より正確には確率であり、
0.8であれば、そのボクセルの80%がおそらく灰白質ですよということになります。
そして、算出の仕方ですが、VBMの場合、
まずは、segmentationを行った際に、灰白質画像、白質画像などのprobability (or density) の計算が先になります。
その後、非線形変換を行う際にjacobian determinantなどが計算され、
その値を用いてmodulationがかけられます。
具体的なところは、拙著すぐできるVBMのmodulationのところをご覧いただけたらと思います。
以上、よろしくお願いします。
この度はご指導有難うございました。bounding boxとvoxel sizeを変更できるスクリプトがありましたので御報告させて頂きます。
http://www-personal.umich.edu/~nichols/JG5/resize_img.m
高橋先生
ありがとうございます。
Ged Ridgwayがそのスクリプトをさらに改良したresize_img.mを発表しています。
ご参考まで。
http://www0.cs.ucl.ac.uk/staff/g.ridgway/vbm/resize_img.m
有難うございます。今後こちらをぜひ使用させて頂きたいと思います。
お忙しいところ、丁寧にご教授頂きまして、ありがとうございます。
また、お返事が大変遅くなり大変失礼致しました。
扁桃体の容積を測定しなくてはならず、大変困っていたところでした。
本当にありがとうございました。
包括脳チュートリアルでは、いつも大変貴重な情報をありがとうございます。
New Segmentの部分について質問させて頂きたく存じます。
New segmentでsmoothing sizeを聞いてくる欄が見当たらないのですが・・。
またNew segmentではき出される画像はsmwc1xxx.niiではなくrc1xxx.niiです。
通常のsegmentではsmwc1xxx.niiが出てきます。
（SPM8でアップデートは最新です。）
ご多忙のこととは存じますが、何卒よろしくお願い致します。
ご指摘ありがとうございました。
私の言葉足らずでした。
> New segmentでsmoothing sizeを聞いてくる欄が見当たらないのですが・・。
> またNew segmentではき出される画像はsmwc1xxx.niiではなくrc1xxx.niiです。
New Segmentの後、DARTELを使ってNormalize to MNI spaceをする時のファイル名ですね。
New Segmentの時は、
Segmentationで、Tissue -> Native -> DARTEL imported
これで、rc1xxx.niiが作成されます。
その後、DARTELを走らせ、
Normalize to MNIで、Preserve Amountとすると、smwc1xxx.niiとなります。
その情報を補足させていただきます。
ご指摘ありがとうございました。
いつも非常にためになる情報をいただき大変ありがとうございます。
この記事で教えていただいたget_totals.mでマスクを指定して測定した結果と、spm_mask.mでマスクをあらかじめした画像をget_totalsで測定した結果が少し異なります。また、左右反対になっているようです。
VBMの結果からは左>右になりそうな領域の結果が、前者では左右となるため、後者の方が正しいのではないかと思うのですが、少し不安になりコメントをさせていただきました。
VBM8で分画化したm0wrp1*.niiやmwrp1*.niiをWFU PickAtlasで作成したROIでマスクしました。
SPMのメーリングリストでget_totalsを検索しても同様の問題点は議論されていませんでしたが、よく知られた問題なのでしょうか？
nktk様
とても重要なご指摘をありがとうございます。
検討しなければいけないことと思いますが、一応、確認ですが、画像はマスク画像も含めてすべてNIfTI形式ですよね？
NIfTIでいじるときに、左右がひっくり返ることはちょっと想像しがたいのですが、検証が必要かと思います。
少し時間をいただいて、私の方でも検証させていただきたいと思います。
よろしくお願いします。
早速のお返事ありがとうございます。
以下、私の環境を少し補足させていただきます。
ファイル形式はすべてNIfTI形式です。
MRIcroやSPMでの表示は、ROI画像の左右は間違っていませんでした。
get_totalsは
http://www0.cs.ucl.ac.uk/staff/g.ridgway/vbm/get_totals.m
のものをそのまま使用しています。
最初のコメントの中で脱字がありましたので、冗長となりますが念のため訂正させていただきます。
誤：VBMの結果からは左>右になりそうな領域の結果が、前者では左右となるため、後者の方が正しいのではないかと思うのですが、少し不安になりコメントをさせていただきました。
正：VBMの結果からは左が右より大きくなりそうな領域の結果が、前者では左が右より小さくなるため、後者の方が正しいのではないかと思うのですが、少し不安になりコメントをさせていただきました。
お忙しいところ、お手間を増やして本当に申し訳ありません。どうぞよろしくお願いします。
nktk様
どうもありがとうございます。
今、状況を再現できました。
確かにマスクがひっくり返ってしまうような問題が起こるようですね。
Ged Ridgwayに質問してみます。
貴重なご指摘をありがとうございました。
nktk様
Ged Ridgwayから返事をいただいた結果、画像の座標系の問題であることがわかりました。
解決策が見つかりましたので、ブログの内容を更新させていただきました。
これで試していただけたらと思います。よろしくお願いします。