【FSL】FSLを用いた画像の位置合わせ ~Registration~


1. 目的
2. 位置合わせで使用する変換について
2.1. 拡大縮小・回転・平行移動・せん断を用いた変換
2.2. 非線形変換
3. コマンド
3.1. FLIRT
3.2. FNIRT
4. 使用例
4.1. 同一被験者における脳画像の位置合わせ(剛体変換)
4.2. 異なる被験者脳の位置合わせ(アフィン変換)
4.3. 異なる被験者脳の位置合わせ(アフィン変換+非線形変換)
4.4. 標準空間上にあるラベル(関心領域)を個人脳に位置合わせ(アフィン変換+非線形変換)


1. 目的

  • 位置合わせで使用する種々の変換の理解
  • 拡大縮小・回転・平行移動・せん断を用いた変換
  • 非線形変換

2. 位置合わせで使用する変換について

位置合わせで使用する変換として、拡大縮小・回転・平行移動・せん断を用いた変換と非線形変換がある。

2.1. 拡大縮小・回転・平行移動・せん断を用いた変換

平行移動および回転を組み合わせた変換を剛体変換(Rigid transform)、拡大縮小・回転・せん断のようにY=AXで表現できる変換を線形変換(1次変換, linear transform)という。また、線形変換に平行移動を組み合わせたY=AX+Bの変換を、アフィン変換(Affine transform)という。詳細は、こちらを参考にすると分かりやすい。

脳MRI画像のような、3次元データの位置合わせの場合、拡大縮小・回転・平行移動・せん断それぞれの自由度は3である。つまり、剛体変換の自由度は6線形変換の自由度は9アフィン変換の自由度は12となる。

2.2. 非線形変換

脳のしわや脳室等の形は個人差があるので、個人脳を標準脳に合わせる際に、拡大縮小・回転・平行移動・せん断を用いた変換では、十分に位置合わせができない。そこで、脳のしわや脳室等までも合わせるために、非線形変換を導入する。

非線形変換は、Y=AX+Bのような線形式で表現できない変換であり、アルゴリズムとしてB-spline法がよく用いられている。詳細は、こちらを参考にするとよい。

3. コマンド

FSLコマンドの、flirtで拡大縮小・回転・平行移動・せん断、fnirtで非線形変換を実行することができる。

3.1. FLIRT

flirtのヘルプは次の通り。

Usage: flirt [options] -in <inputvol> -ref <refvol> -out <outputvol>
       flirt [options] -in <inputvol> -ref <refvol> -omat <outputmatrix>
       flirt [options] -in <inputvol> -ref <refvol> -applyxfm -init <matrix> -out <outputvol>

  Available options are:
        -in  <inputvol>                    (no default)
        -ref <refvol>                      (no default)
        -init <matrix-filname>             (input 4x4 affine matrix)
        -omat <matrix-filename>            (output in 4x4 ascii format)
        -out, -o <outputvol>               (default is none)
        -datatype {char,short,int,float,double}                    (force output data type)
        -cost {mutualinfo,corratio,normcorr,normmi,leastsq,labeldiff,bbr}        (default is corratio)
        -searchcost {mutualinfo,corratio,normcorr,normmi,leastsq,labeldiff,bbr}  (default is corratio)
        -usesqform                         (initialise using appropriate sform or qform)
        -displayinit                       (display initial matrix)
        -anglerep {quaternion,euler}       (default is euler)
        -interp {trilinear,nearestneighbour,sinc,spline}  (final interpolation: def - trilinear)
        -sincwidth <full-width in voxels>  (default is 7)
        -sincwindow {rectangular,hanning,blackman}
        -bins <number of histogram bins>   (default is 256)
        -dof  <number of transform dofs>   (default is 12)
        -noresample                        (do not change input sampling)
        -forcescaling                      (force rescaling even for low-res images)
        -minsampling <vox_dim>             (set minimum voxel dimension for sampling (in mm))
        -applyxfm                          (applies transform (no optimisation) - requires -init)
        -applyisoxfm <scale>               (as applyxfm but forces isotropic resampling)
        -paddingsize <number of voxels>    (for applyxfm: interpolates outside image by size)
        -searchrx <min_angle> <max_angle>  (angles in degrees: default is -90 90)
        -searchry <min_angle> <max_angle>  (angles in degrees: default is -90 90)
        -searchrz <min_angle> <max_angle>  (angles in degrees: default is -90 90)
        -nosearch                          (sets all angular search ranges to 0 0)
        -coarsesearch <delta_angle>        (angle in degrees: default is 60)
        -finesearch <delta_angle>          (angle in degrees: default is 18)
        -schedule <schedule-file>          (replaces default schedule)
        -refweight <volume>                (use weights for reference volume)
        -inweight <volume>                 (use weights for input volume)
        -wmseg <volume>                    (white matter segmentation volume needed by BBR cost function)
        -wmcoords <text matrix>            (white matter boundary coordinates for BBR cost function)
        -wmnorms <text matrix>             (white matter boundary normals for BBR cost function)
        -fieldmap <volume>                 (fieldmap image in rads/s - must be already registered to the reference image)
        -fieldmapmask <volume>             (mask for fieldmap image)
        -pedir <index>                     (phase encode direction of EPI - 1/2/3=x/y/z & -1/-2/-3=-x/-y/-z)
        -echospacing <value>               (value of EPI echo spacing - units of seconds)
        -bbrtype <value>                   (type of bbr cost function: signed [default], global_abs, local_abs)
        -bbrslope <value>                  (value of bbr slope)
        -setbackground <value>             (use specified background value for points outside FOV)
        -noclamp                           (do not use intensity clamping)
        -noresampblur                      (do not use blurring on downsampling)
        -2D                                (use 2D rigid body mode - ignores dof)
        -verbose <num>                     (0 is least and default)
        -v                                 (same as -verbose 1)
        -i                                 (pauses at each stage: default is off)
        -version                           (prints version number)
        -help

基本的な使い方は、以下。

単純に、位置合わせを実行したい場合。

flirt -in <位置合わせしたい画像> -ref <位置合わせ先の画像> -dof <自由度> -out <出力画像>

一度、変換行列を生成して、次にそれを適応する場合。

# 変換行列を生成
flirt [options] -in <位置合わせしたい画像> -ref <位置合わせ先の画像> -dof <自由度> -omat <変換行列の出力ファイル>

# 変換行列を適用
flirt [options] -in <位置合わせしたい画像> -ref <位置合わせ先の画像> -applyxfm -init <適用したい変換行列> -out <出力画像>

3.2. FNIRT

fnirtのヘルプは次の通り。

Usage: 
fnirt --ref=<some template> --in=<some image>
fnirt --ref=<some template> --in=<some image> --infwhm=8,4,2 --subsamp=4,2,1 --warpres=8,8,8

Compulsory arguments (You MUST set one or more of):
	--ref		name of reference image
	--in		name of input image

Optional arguments (You may optionally specify one or more of):
	--aff		name of file containing affine transform
	--inwarp	name of file containing initial non-linear warps
	--intin		name of file/files containing initial intensity mapping
	--cout		name of output file with field coefficients
	--iout		name of output image
	--fout		name of output file with field
	--jout		name of file for writing out the Jacobian of the field (for diagnostic or VBM purposes)
	--refout	name of file for writing out intensity modulated --ref (for diagnostic purposes)
	--intout	name of files for writing information pertaining to intensity mapping
	--logout	Name of log-file
	--config	Name of config file specifying command line arguments
	--refmask	name of file with mask in reference space
	--inmask	name of file with mask in input image space
	--applyrefmask	Use specified refmask if set, default 1 (true)
	--applyinmask	Use specified inmask if set, default 1 (true)
	--imprefm	If =1, use implicit masking based on value in --ref image. Default =1
	--impinm	If =1, use implicit masking based on value in --in image, Default =1
	--imprefval	Value to mask out in --ref image. Default =0.0
	--impinval	Value to mask out in --in image. Default =0.0
	--minmet	non-linear minimisation method [lm | scg] (Levenberg-Marquardt or Scaled Conjugate Gradient)
	--miter		Max # of non-linear iterations, default 5,5,5,5
	--subsamp	sub-sampling scheme, default 4,2,1,1
	--warpres	(approximate) resolution (in mm) of warp basis in x-, y- and z-direction, default 10,10,10
	--splineorder	Order of spline, 2->Quadratic spline, 3->Cubic spline. Default=3
	--infwhm	FWHM (in mm) of gaussian smoothing kernel for input volume, default 6,4,2,2
	--reffwhm	FWHM (in mm) of gaussian smoothing kernel for ref volume, default 4,2,0,0
	--regmod	Model for regularisation of warp-field [membrane_energy bending_energy], default bending_energy
	--lambda	Weight of regularisation, default depending on --ssqlambda and --regmod switches. See user documentation.
	--ssqlambda	If set (=1), lambda is weighted by current ssq, default 1
	--jacrange	Allowed range of Jacobian determinants, default 0.01,100.0
	--refderiv	If =1, ref image is used to calculate derivatives. Default =0
	--intmod	Model for intensity-mapping [none global_linear global_non_linear local_linear global_non_linear_with_bias local_non_linear]
	--intorder	Order of polynomial for mapping intensities, default 5
	--biasres	Resolution (in mm) of bias-field modelling local intensities, default 50,50,50
	--biaslambda	Weight of regularisation for bias-field, default 10000
	--estint	Estimate intensity-mapping if set, default 1 (true)
	--numprec	Precision for representing Hessian, double or float. Default double
	--interp	Image interpolation model, linear or spline. Default linear
	-v,--verbose	Print diagnostic information while running
	-h,--help	display help info

基本的な使い方は、以下。

単純に、位置合わせしたい場合。

fnirt --ref=<位置合わせ先の画像> --in=<位置合わせしたい画像> --iout=<出力画像>

最初にflirtを実行し、次にfnirtを実行して、位置合わせする場合。

flirt [options] -in <位置合わせしたい画像> -ref <位置合わせ先の画像> -dof <自由度> -omat <flirt変換行列の出力ファイル>
fnirt --in=<位置合わせしたい画像> --aff=<適用したいflirt変換行列>  --cout=<fnirt変換行列の出力ファイル>

applywarp --in=<位置合わせしたい画像> --ref=<位置合わせ先の画像> --warp=<fnirt変換行列の出力ファイル> --out=<flirt/fnirt後の出力画像>

fnirtでは、様々なパラメータを設定することができるが、FSLでは種々のパラメータ値が記載された設定ファイル(.cnf)を提供している。設定ファイルは、${FSLDIR}/etc/flirtschで確認することができる。fnirtで設定ファイルを用いるには、-configオプションを用いて設定ファイルを指定する。

${FSLDIR}/etc/flirtsch
├── FA_2_FMRIB58_1mm.cnf  # 個人FAとFMRIB58_1mm(標準FA)の設定ファイル
├── GM_2_MNI152GM_2mm.cnf  # 個人灰白質(GM)とMNI152GM_2mm(標準GM)の設定ファイル
├── T1_2_MNI152_2mm.cnf  # 個人T1WIとMNI152_2mm.cnf(標準T1WI)の設定ファイル
├── b02b0.cnf  # 拡散強調画像のb=0画像の位置合わせで用いる設定ファイル
├── b02b0_1.cnf  # 拡散強調画像のb=0画像の位置合わせで用いる設定ファイル
├── b02b0_2.cnf  # 拡散強調画像のb=0画像の位置合わせで用いる設定ファイル
└── b02b0_4.cnf  # 拡散強調画像のb=0画像の位置合わせで用いる設定ファイル

4. 使用例

flirtおよびfnirtを用いた脳画像の位置合わせについて、使用例を用いて解説していく。

位置合わせの精度を高めるために、位置合わせの前に前処理として頭蓋除去をしておくとよい。やり方は、以下の記事を参考にするとよい。

頭蓋除去をすると次のようになる。

4.1. 同一被験者における脳画像の位置合わせ(剛体変換)

標準空間上にある個人脳(T1_skull_stripped_inMNI.nii.gz)を、個人空間上の個人脳(T1_skull_stripped.nii.gz)に位置合わせする。

二つの画像には、個人脳と標準脳のとの間には、この程度の位置ずれがある。

標準空間上にある個人脳(T1_skull_stripped_inMNI.nii.gz)を、個人空間上の個人脳(T1_skull_stripped.nii.gz)に位置合わせするには、次のコマンドを実行する。同一被験者脳の位置合わせであるため、自由度(DoF)は6としている。

flirt -in T1_skull_stripped_inMNI.nii.gz -ref T1_skull_stripped.nii.gz -dof 6 -out T1_skull_stripped_MNI2individual.nii.gz

標準空間から個人空間に位置合わせした個人脳と、個人空間にある個人脳を重ね合わせると、次のようになる。

4.2. 異なる被験者脳の位置合わせ(アフィン変換)

頭蓋除去済みのT1WI(T1_skull_stripped.nii.gz)を標準脳T1WI(MNI152_T1_1mm_brain.nii.gz)に位置合わせする。MNI152_T1_1mm_brain.nii.gzは、「${FSLDIR}/data/standard/MNI152_T1_1mm_brain.nii.gz」にある。

個人脳と標準脳のとの間には、この程度の位置ずれがある。

個人脳を標準脳に合わせるには、次のコマンドを実行する。この時、自由度は12としている。

flirt -in T1_skull_stripped.nii.gz -ref MNI152_T1_1mm_brain.nii.gz -dof 12 -out T1_skull_stripped_inMNI.nii.gz

標準脳に位置合わせした個人脳と標準脳を重ね合わせると、次のようになる。

4.3. 異なる被験者脳の位置合わせ(アフィン変換+非線形変換)

頭蓋除去済みのT1WI(T1_skull_stripped.nii.gz)を標準脳T1WI(MNI152_T1_1mm_brain.nii.gz)に位置合わせする場合、アフィン変換だけでは、脳のしわや脳室等における位置合わせ不十分である(下図)。

そこで、アフィン変換にあわせて非線形変換も組み合わせる。

アフィン変換および非線形変換を用いて、頭蓋除去済みのT1WI(T1_skull_stripped.nii.gz)を標準脳T1WI(MNI152_T1_1mm_brain.nii.gz)に位置合わせするには、以下のコマンドを実行する。

flirt -in T1_skull_stripped.nii.gz -ref MNI152_T1_1mm_brain.nii.gz -omat indiv2std.mat
fnirt --ref=MNI152_T1_1mm_brain.nii.gz --in=T1_skull_stripped.nii.gz --aff=indiv2std.mat --iout=T1_skull_stripped_inMNI.nii.gz

アフィン変換および非線形変換を用いた、位置合わせの結果は以下。

4.4. 標準空間上にあるラベル(関心領域)を個人脳に位置合わせ(アフィン変換+非線形変換)

何らかのアトラスで定義された関心領域(ROI)を用いて、個人脳の何らかの定量値を計測したい場合がある。その時には、アトラスを個人脳に位置合わせしなくてはならない。

標準脳のFAにあるアトラス(JHU-ICBM-labels-1mm.nii.gz)を、個人脳のFAに位置合わせする。JHU-ICBM-labels-1mm.nii.gzは、「${FSLDIR}/data/atlases/JHU/JHU-ICBM-labels-1mm.nii.gz 」にある。

手順は、次の通り。

  1. flirtを用いて個人脳FAを標準脳FAに位置合わせ
  2. 1の変換行列を初期値として、fnirtで個人脳FAを標準脳FAに位置合わせ
  3. 2で得られた変換行列(個人脳→標準脳)を反転させて、標準脳から個人脳へと変換する行列に
  4. 標準脳上にあるアトラスに、3の変換行列を適用して、アトラスを個人脳に位置合わせ
# 位置合わせ
flirt -in FA.nii.gz -ref FMRIB58_FA_1mm.nii.gz -dof 12 -omat indiv2std.mat
fnirt --in=FA.nii.gz --aff=indiv2std.mat --config=FA_2_FMRIB58_1mm.cnf --cout=warp_indiv2std.nii.gz

# 変換行列(個人脳→標準脳)を反転して、標準脳から個人脳へと変換する行列に
invwarp -w warp_indiv2std.nii.gz -o warp_std2indiv.nii.gz -r FA.nii.gz

# 変換行列を適用して、アトラスを個人脳に位置合わせ
applywarp --in=JHU-ICBM-labels-1mm.nii.gz --ref=FA.nii.gz --warp=warp_std2indiv.nii.gz --interp=nn --out=JHU-ICBM-labels-1mm_indiv.nii.gz

位置合わせの結果は、次の通り。

Print Friendly, PDF & Email

【FSL】FSLを用いた画像の位置合わせ ~Registration~” へのコメント

  1. Lin4neuro、いつもお世話になっております。
    唐突ではございますが、同一患者の3D FLAIR とT2*強調像を重ねた画像を作成することは可能でしょうか。

    • 可能です。

      以下のように行います。

      – 3D FLAIR に対して bet で skull-stripping を行います。入力ファイルが flair.niiとし、flair_brain.nii.gz という出力ファイル名にするならば以下のようになります。

      bet flair.nii flair_brain.nii.gz
      

      その次に、t2* と flair_brain を FLAIR で位置合わせします。剛体変換なので、-dof 6 とするのがいいでしょう。t2* の入力を t2star.nii, 位置合わせをしたものをt2star_realigned.nii.gz とします。

      flirt -in t2star.nii -ref flair_brain.nii.gz -dof 6 -out t2star_realigned.nii.gz
      

      これでできると思いますがいかがでしょうか。

      • ありがとうございます。
        操作自体は無事にできたようですが、T2*強調像が目立っていてCSFが高信号です。
        参考文献の画像ではFLAIR寄りの画像でしたので、雰囲気が異なります。
        別の方法もありましたら、ご教示いただけると幸いです。

          • ご返信ありがとうございます。
            t2star_realigned.nii.gzをFSLeyesで見ると、T2*強調像のみしか映っていないように見えます。重なっているかどうか判断しかねます。
            もっと大脳白質病変がFLAIRと同様に高信号に見えると良いのですが…

            ちなみに、おしえて頂いた、
            『bet flair.nii flair_brain.nii.gz』
            flair.nii は flair.nii.gz ではないということでよろしかったでしょうか。
            『flirt -in t2star.nii -ref flair_brain.nii.gz -dof 6 -out t2star_realigned.nii.gz

            t2star.nii は t2star.nii.gz ではないということでよろしかったでしょうか。

          • NIFTI画像は nii 形式と nii.gz 形式があります。niiをgz形式で圧縮したものが nii.gz 形式です。
            FSLは nii, nii.gz どちらも扱えるので、どちらでも構わないです。

          • 承知しました。
            ご対応いただきありがとうございます。

  2.  根本先生

     いつもお世話になっています。またひとつ、質問させてください。
     
     今回の記事の最後にinvwarpを用いて、変換を反転させていらっしゃいます。FLIRT+FNIRTの場合にはFNIRTのcoutで作成したNIFTI画像を使っていらっしゃいますが、FLIRT単独の変換で終わる場合にはどのようにすればよいのでしょうか…? 
     お手すきの際にご教示いただけましたら幸いです。

    • 赤池先生

      FLIRT単独で何をしたいと考えていらっしゃいますか?

      それによって変わるかなと思います。

      •  根本先生

         早々にありがとうございます。
         MNI空間で作成したマスクを各症例のFLAIR像にあわせたいというのが最終的な目的です。

         多発性硬化症などの病変がある脳画像のFLAIR像から、各部位の病変量(前頭葉には●mlなどの形で)を算出したいと考えています。SPMのtoolboxを用いて、各症例のFLAIR像から病変マスクを造ることができましたので、MNI空間のラベルからmaskを作成し、各症例のFLAIR像にあわせて変換し、病変マスクと合成しようと考えました。
         FLAIR像を標準脳に変換するのに、個人脳FLAIR→個人脳T1WI→標準脳に変換という方法をとりました。標準脳→個人脳FLAIRの行列を作成したかったのですがうまくいきませんでした。ブログの記事を参考に標準脳を個人脳T1WIに変関する方法はわかりましたので、そこからFLAIRにあわせる(つまり、FLIRTの変換行列にinvwarpを適用する)ことを試みたのですが、こちらもうまく変換できず、質問させていただいたという次第です。
         

        • そういうことですね。

          それならば、

          flirt -dof 6 -in FLAIR画像 -ref T1画像 -out FLAIR画像_r

          などとすると、T1画像に位置合わせをしたFLAIR画像が得られます。ファイル名の後ろに _r がついています。

          それは、T1と重ね合わせがなされていますので、

          maskを invwarp すれば、FLAIR画像とあうはずです。

          それでいかがでしょうか?

          •  いわれてみましたら、仰る通りの方法でうまくいくはずですね…すみません、お騒がせいたしました。
             ただ、できれば、すでに作成した病変マスク(元画像のFLAIRから作成)を使用したいのですが、FLIRTで得られた変換行列をapplywarpで使用することは可能なのでしょうか?

          • 以下でどうでしょう。

            flirt -in FLAIR画像 -ref T1画像 -out FLAIR画像_r -omat FLAIR2T1.mat -dof 6
            flirt -interp nearestneighbour -in mask画像 -ref T1画像 -out mask画像_r -init FLAIR2T1.mat -applyxfm

          •  根本先生

             ありがとうございます。
             このような使い方があるのですね…勉強になります。ぜひ、試してみたいと思います。

          •  根本先生

             教えていただいた方法でうまくいきました!
             ありがとうございます。

  3. ピングバック: 【FSL】Voxel-BasedAnalysis: VBA

コメントを残す

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