VBMの結果を用いて容積を測定する方法

VBMで得られる画像を用いて、灰白質や白質画像の容積を求めることができます。

事前の準備として、以下のことが必要です。

  • get_totals.m
  • Ged Ridgwayが公開している容積算出のためのスクリプトです。http://www0.cs.ucl.ac.uk/staff/g.ridgway/vbm/get_totals.mここから入手できます(右クリックしてリンク先を保存としてください)。これをSPMのディレクトリに保存します。

  • WFU PickAtlas
  • これは、局所容積を求めたい時に便利です。マスク画像を簡単に作成することができます。http://fmri.wfubmc.edu/software/PickAtlasここから入手できます。今回はこのセットアップ方法は説明しませんのであしからず。

画像としては、以下のものが必要となります。

  • SPM8のNew Segmentを使っている時
  • modulationをかけ、smoothを[0 0 0]としたもの。出力ファイルは、灰白質ならば、smwc1xxx.niiのような形をとるはずです。New Segmentでは、smoothingをかけないという選択肢はないので、[0 0 0]とします。今までいろいろ試してみましたが、smoothingをかけると信号値が若干変化しますので、smoothing前のデータを用いるのがよいと思います。

  • SPM12のNew Segmentを使っている時
  • modulationをかけ、smoothを[0 0 0]としたもの。出力ファイルは、灰白質ならば、mwc1xxx.niiのような形をとるはずです。SPM8との違いは、smoothingを[0 0 0]とすると、sがつかないようになったようです。(猪狩先生、情報をどうもありがとうございます)

  • VBM8を使っている時
  • Modulationを”affine + non-linear (SPM8 default)”とします。non-linear onlyですと、線形変換の結果は考慮されませんので、頭蓋容積を補正済みのデータになってしまい、生値を得ることができないので注意が必要です。出力ファイルは、灰白質ならば、mwrp1xxx.niiとなります。

それでは、はじめましょう。

  1. 基本編
  2. まずは、全灰白質容積を算出しましょう。

    Matlabのコマンドウィンドウから、

    get_totals
     
    

    とタイプします。そうすると、画像を選択するダイアログが出ますので、灰白質容積を求めたい画像を選択します。(複数選択できます。)

    すると、Matlabのコマンドウィンドウに、

    ans =
      739.8652
    

    というように出力されます。これが求めたい容積です。複数画像を選択した場合、ファイルの選択順に出力されます。

    これは簡単ですね。get_totalsとコマンドを打てばいいのですから。

    それでは、次にいきましょう。

  3. 応用編
  4. 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)”を選択します。

      coregister_reslice

      そうすると、下図のようなウィンドウが出現しますので、

      Image Defining Spaceに、灰白質画像を指定します。
      Images to Resliceに、マスク画像を指定します。
      InterpolationはNearest neighbourにします。そうすることで、マスクがバイナリ画像(二値画像)のままとなります。

      その他の設定はデフォルトのままで結構です。

      reslice_settings

      そうすると、計算が始まり、新しくファイルの先頭にrが付加されたファイルが生成されます。今の場合ですと、rlt_hemisphere.niiという画像ができます。
      これをSPMのDisplayで表示してみます。そうすうと、Vox Sizeが-1.5 x 1.5 x 1.5と変わっていることがわかるかと思います。

      ここまで来ればあと一息です。

    • get_totalsで左半球の灰白質容積算出
    • Matlabのコマンドウィンドウで、次のようにタイプします。

      get_totals([],[],[])
       
      

      ここで、[]はMatlabではそこが空ということを意味します。今は、画像ファイルも閾値もマスク画像も空ですよという意味です。そうすると、最初に下図のように画像ファイルを選択するダイアログが表示されます。

      そうすると、次にマスク画像を選択する画面が出てきます。ここでは、先程準備したrlt_hemisphere.niiを選択します。

      なお、閾値は空ですと、0に設定されます。もし、信号値0.1以上を示す領域での容積を求めたいならば、ここで[]の代わりに、0.1と入力すれば大丈夫です。

      そうすると、Matlab上に

      ans =
        350.6627
      

      のような表示が出てきます。先ほどの全灰白質容積と比較すると約半分になっていることがわかりますね。

このようにしてVBMで得られた画像から容積を求めることができます。

Print Friendly, PDF & Email

27 thoughts on “VBMの結果を用いて容積を測定する方法

  1. いつも本サイトや根本先生の著書「すぐできる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,[])で大丈夫です。マスク画像については、ご自身で決める必要があります。通常は脳をカバーするマスク画像を使うことが多いかなと思います。

      以上、よろしくお願いします。

      根本清貴

      • 根本先生

        ご返信ありがとうございます.2に関して,ボクセル値=定量値の場合,get_totalsから求まる値は,volumeではなくて”定量値の合計”という理解でよろしいでしょうか.この場合,”volume”や”定量値の平均”を求める方法はございますでしょうか.

        打田佑人

        • 打田先生

          ボクセル値=定量値の場合,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の値を出力してくれます)

            定量画像などの場合は、それぞれ慎重に検討する必要があると思います。

            根本清貴

          • 根本先生

            ありがとうございます.最後にすみません,先生は,「体積」と「容積」をどのように使い分けておられますでしょうか.

            打田佑人

          • 打田先生

            すみません、使い分けていません…。
            その場で頭に浮かんできた方を使っているだけなので気にされないでください。

            根本清貴

          • 根本先生

            細かい質問になり申し訳ありません.根本先生の説明はいつも大変わかりやすくて助けられています.今後ともどうぞよろしくお願いいたします.

            打田佑人

  2. いつも先生のサイトを読ませて頂き、研究を行なっております。
    東京医科歯科大 耳鼻科 山本と申します。
    VBMで得られた画像より、get_totallsを使用したりしてROIの体積などは解析しております。
    VBMでは密度も算出できるとなっておりますが、具体的にどの数値になるのかご教授頂けますでしょうか?
    得られた体積を使用して、算出するのでしょうか?
    基本的な質問で申し訳ないことでございます。

    • 山本先生

      すみません、先生のコメントがなぜかスパム扱いされており、今しがた気づきました。

      VBMでの密度(density)は、純粋に0-1の間となります。より正確には確率であり、
      0.8であれば、そのボクセルの80%がおそらく灰白質ですよということになります。

      そして、算出の仕方ですが、VBMの場合、

      まずは、segmentationを行った際に、灰白質画像、白質画像などのprobability (or density) の計算が先になります。
      その後、非線形変換を行う際にjacobian determinantなどが計算され、
      その値を用いてmodulationがかけられます。
      具体的なところは、拙著すぐできるVBMのmodulationのところをご覧いただけたらと思います。

      以上、よろしくお願いします。

  3. お忙しいところ、丁寧にご教授頂きまして、ありがとうございます。
    また、お返事が大変遅くなり大変失礼致しました。

    扁桃体の容積を測定しなくてはならず、大変困っていたところでした。
    本当にありがとうございました。

  4. 包括脳チュートリアルでは、いつも大変貴重な情報をありがとうございます。

    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となります。

      その情報を補足させていただきます。

      ご指摘ありがとうございました。

  5. いつも非常にためになる情報をいただき大変ありがとうございます。
    この記事で教えていただいた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から返事をいただいた結果、画像の座標系の問題であることがわかりました。
            解決策が見つかりましたので、ブログの内容を更新させていただきました。
            これで試していただけたらと思います。よろしくお願いします。

打田佑人 へ返信する コメントをキャンセル

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください