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

【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 」にある。

手順は、次の通り。

flirtを用いて個人脳FAを標準脳FAに位置合わせ
1の変換行列を初期値として、fnirtで個人脳FAを標準脳FAに位置合わせ
2で得られた変換行列（個人脳→標準脳）を反転させて、標準脳から個人脳へと変換する行列に
標準脳上にあるアトラスに、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

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

著者情報：斎藤勇哉

順天堂大学医学部大学院医学研究科放射線診断学講座所属
脳MRI 画像解析が専門であり、テーマは①神経変性疾患の機序解明、②医用人工知能の開発、③多施設データのハーモナイゼーション、④速読が脳に与える影響や学習効果、⑤SNS解析を用いたマーケティング戦略の改善。
医療分野に関わらず、自然言語処理・スクレイピング・データ分析・Web アプリ開発を得意とし、企業や他大学の研究を支援。
主な使用言語は、Python、Shell Script、MATLAB、HTML、CSS