FSLを用いたシェルスクリプト演習:脳3次元T1強調MRIの水平断一覧画像の作成

あるプロジェクトにおいて、「MRI画像の一覧画像を作成してもらいたい」という依頼を受けました。

イメージとしては、下記のようなものです。

LB_example

以前、FSLのシェルスクリプトチュートリアルを翻訳していた際に、演習に
Lightbox Viewerがあったのを思い出して、改めて、それを見ながら、さらに機能を追加したスクリプトを考えてみました。

作成しながら勉強になることがいくつもあったので、細かくわけてみたいと思います。

  1. 剛体変換 (Rigid body transformation)
  2. MRIは、人によって頭の位置が異なりますので、ある範囲内で画像を表示しようとすると、脳が変なところで切れてしまうことがあります。これを解決するために、個々人の脳を、標準脳に対して線形変換のうちの剛体変換 (rigid body transformation)を行おうと考えました。剛体変換は、脳を水平移動と回転だけするものです。つまり、脳の大まかな形は変わりません。水平移動と回転をx, y, z軸の3方向に対して行いますので、2(水平移動、回転)×3軸で、6つのパラメーターを使うことになります。

    FSLでは、flirtというプログラムが線形変換を行うプログラムです。用法としては、以下になります。

    flirt [オプション] -in 入力画像 -ref 参照画像 -out 出力画像
    

    今、オプションに、剛体変換を指定したいと思います。flirtでは、-dof の後に数値を指定します。今は上記のように6パラメーターですので、-dof 6となります。

    さらに、参照画像ですが、FSLの中に入っているMNI152_T1_1mm.nii.gzを使いたいと考えました。このファイルは、FSLの中のdata/standard/ディレクトリに入っています。FSLのパスは、${FSLDIR}で指定できます。これを使うことで、MacでもLinuxでも使えるスクリプトになります。

    あとで、スクリプトを読みやすくするために、参照画像は変数 $ref で指定することにします。そして、出力画像は入力画像のファイル名の最初に”r”を付加することにします。

    入力画像はforループを使って変数fにおさめられているとします。

    すると、以下のようになります。

    ref=${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz
    
    for f in "$@"
    do
      #Rigid body registration
      flirt -dof 6 -in $f -ref $ref -out r${f}
    (以下続く)
    

    これで、MNI152_T1_1mm.nii.gzと同じFOVの画像を得ることができます。

    具体的に見てみましょう。

    元々はこんな画像です。

    orig_t1

    これが、出力画像として、こんな画像になります。

    flirted_t1

    首の部分がきれいに切れていることがわかります。また、AC-PCラインに沿って脳の位置がpitch方向に回転していることがわかります。

  3. 水平断の抽出
  4. 次に線形変換後の画像から水平断の画像を抽出しようと考えました。線形変換後画像を見ると、z=26からz=145に脳がだいたい入っていることがわかりました。そして、5×6で表示したいと思っていましたので、30枚抽出できたらいいなと考えました。z=26からz=145だと120枚になりますので、4枚おきにすれば30枚選ぶことができます。

    MRI画像から水平断を抽出し、PNG画像として保存するコマンドにFSLにはslicerというコマンドがあります。(3D-Slicerとは違うので注意が必要です)

    水平断を抽出する用法としては、以下になります。

    slicer 入力画像 -z -番号 出力画像
    

    今、これを、z=26から4枚おきに画像を生成したいと思います。

    4つごとに数字を出す方法に、UNIX系OSでは、seqというコマンドがあります。

    seq 開始番号 間隔値 終了番号
    

    なお、seqには-fオプションを使うことで、ゼロ埋めを行うことができます。つまり、26は026と表すことができます。これは、”%03g”となるそうです。

    つまり、

    seq -f "%03g" 26 4 145
    

    となります。

    今、出力画像は、ファイル名の後に連番をつけたいと思います。
    ただ、ファイル名の拡張子はとりたいと思います。FSLには、remove_extというプログラムがあります。

    これらをまとめて、以下のようにすることで、z=26から4枚ごとに30枚分の水平断画像を抽出することができます。バッククオートを使っていることに注意してください。

    for i in `seq -f "%03g" 26 4 145`
      do
       file=`remove_ext $f`
       output=${file}_$i.png
       slicer r$f -z -$i $output
      done
    

    これによって30枚のPNG画像を得ることができます。

  5. PNG画像を配列し、1枚の画像とする
  6. 次に、この得られた30枚のPNG画像を5×6のタイル状に配置したいと思います。

    FSLには、pngappendというプログラムがあります。

    pngappendは次のような用法で使います。

    pngappend 入力画像1 +/- 入力画像2 +/- 入力画像3 ... 出力画像
    

    ここで、入力画像を”+”か”-“でつなぐのですが、”+”の時は、水平(横)方向に、”-“の時は垂直(縦)方向につなげることができます。

    つまり、6つの画像(im1〜im6)を横3×縦2としたい場合、以下になります。

    pngappend im1 + im2 + im3 - im4 + im5 + im6 出力画像
    

    今、5×6の配置にしたいので、
    im1 + im2 + im3 + im4 + im5 – im6 + im7 + …

    としないといけません。
    ただ、これを手で書くのは賢くないので、なんとか自動で書くように考えてみます。

    いろいろ試行錯誤したあげく、今回は、以下の方法でいきたいと考えました。

    まず、先ほど作成したPNG画像の一覧を ls で表示します。
    その次に、5行目ごとに、”-” を行頭につけ、その他には、”+”を行頭につけます。
    そして、その後、それらを改行を削除して1つに表示します。
    これらを”inputarg”という変数に代入します。

    ポイントは、「5行目ごとに」となります。

    これにAWKを使うことにしました。
    AWKでは、いくつか決められた変数があり、その中のひとつにNRというものがあります。これは、Number of Recordsの略で、行数を意味します。
    今、6行目、11行目、16行目…にそれぞれ “-” をつけたいと思っています。
    このルールがどこにあるかというと、「5で割ったときに余りが1」となります。

    AWKにおいて、あまりは%で表示されます。今の場合だと、

    NR%5==1
    

    であらわすことができます。

    また、AWKでは、1行1行処理していくわけですが、1行全部は”$0″で表示します。文字列を表示するには、print “文字列”となります。これを合体させると、もし、行数の前に”-“をつけたい場合には、

    print " - " $0
    

    で表示できます。

    これらをまとめて、以下で表示できることがわかりました。

    inputarg=`ls ${file}_???.png | \
    awk '{ if (NR%5==1) {print " - " $0} else {print " + " $0} }'` 
    

    こうすると無事に5行ごとに “-” が、それ以外は “+” がつきます。

    その次に、改行をとりたいと思います。これは、UNIX系OSにある “tr” が便利です。trで改行を削除するには以下を使います。

    tr -d '\n'
    

    そして、最後に一番最初の “-” だけ取り除く必要があります。これは、sedを使うことで簡単にできます。

    sed 's/-//'
    

    これらをまとめると、以下のようになります。

    inputarg=`ls ${file}_???.png | \
    awk '{ if (NR%5==1) {print " - " $0} else {print " + " $0} }' | \
    tr -d '\n' | sed 's/-//'`
    

    出力画像は、ファイル名の前にLB (Lightboxの略とします)をつけることにしました。
    従って、以下になります。

    pngappend $inputarg LB_${file}.png
    
  7. テンポラリファイルの削除
  8. 途中で作成した画像を削除します。

    rm ${file}_???.png
    

以上で、完成です。
若干機能を加えてありますが、以下のようになりました。

#!/bin/bash
#Generate Lightbox View of 3D T1 image

FSLOUTPUTTYPE=NIFTI
ref=${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz

#Check if arguments are specified.
if [ $# -lt 1 ] ; then
 echo "Usage: $0 <filename> [options]" ;
 exit 1 ;
fi

for f in "$@"
do
  #Rigid body transformation
  if [ ! -e r${f} ]; then
    flirt -dof 6 -in $f -ref $ref -out r${f}
  fi

  #Obtain every 4 slices from z=26 to 145
  for i in `seq -f "%03g" 26 4 145`
  do
    file=`remove_ext $f`
    output=${file}_${i}.png
    slicer r$f -z -${i} $output
  done
  
  #Prepare argument for pngappend
  inputarg=`ls ${file}_???.png | \
    awk '{ if (NR%5==1) {print " - " $0} else {print " + " $0} }' | \
    tr -d '\n' | sed 's/-//'`

  #For debug
  #  echo $inputarg

  #Append png files 
  pngappend $inputarg LB_${file}.png
  
  #Delete temporary files
  rm ${file}_???.png
  #rm r$f
done

試してみたい方は、以下のスクリプトをダウンロードして試してみてください。

lightbox-3DT1.shをダウンロード(右クリックで名前をつけて保存で保存してください)

Print Friendly, PDF & Email

4 thoughts on “FSLを用いたシェルスクリプト演習:脳3次元T1強調MRIの水平断一覧画像の作成

  1. 根本先生

    いつも先生のサイトではいろいろと勉強させていただきありがとうございます.
    大学病院で神経内科医をしております.

    質問する場所がここでよいかどうかわからないのですが,本日はFSL-FIRSTを用いたvertex analysisの件でお伺いしたいことがあります.

    FSL-FIRSTでsegmentationした結果を,3D slicerを用いて確認してみると,海馬と扁桃体が接する部分など,いくつかの部位でsegmentationがうまくいっていないことがしばしば見られます.そこで3D slicerを用いて,segmentationがうまくいっていない部位をmanualで修正して,海馬や扁桃体などの体積を算出しております.

    先生にお伺いしたいのは,manual修正した結果を反映させて,vertex analysisを行うことはできるのでしょうか?

    FSLWikiのマニュアルを参照しますと,FSL-FIRSTでsegmentationした際に作成されるbvarsファイルなどを用いてvertex analysisを行うことになっています.しかしそもそも最初のsegmentationがうまくいっていない場合,正確なvertex analysisが行えないのではないかと考えます.

    なにかご存じでしたら教えていただけますでしょうか?

    • ご質問ありがとうございます。

      残念ながら、私はFSL-FIRSTをそこまで使いこなしてはいません。
      岩手医大の山下先生が詳しいので、彼に聞いてみて、その結果をお伝えするようにします。
      よろしくお願いします。

    • 黒川先生

      岩手医大の山下先生に問い合わせたところ、以下の返事をいただきました。

      run_first_allの中身除いてみましたが、メインの計算はrun_firstを呼び出していて、それがさらにfirstという実行ファイルを呼び出しているんですけど、一連の計算が全てこのfirstの中で行われているようで、融通がきかないみたいです。

      引数も見てみましたがダメそうでした・・。あとは誰かがメーリングリストとかQ&Aで質問しているかも知れませんので、それを調べるしか(または本家に対応を求めるしか)ないような気がします。

      ということで、簡単ではなさそうです。
      お役に立てず恐縮ですが、以上、よろしくお願いします。

      • 根本先生

        山下先生にお問い合わせいただき大変ありがとうございました.
        マニュアル修正の結果をvertex analysisに反映するのは現時点では難しそうですね.

        過去のメーリングリストやQ&Aの中身も検索してみようと思います.
        なにかわかりましたらご報告させていただきます.

kiyotaka へ返信するコメントをキャンセル

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