【MRtrix】MRtrixを用いたギブズのリンギングアーチファクト(Gibbs ringing)の除去


1. 目的
2. コマンド
3. 使用例
3.1. 構造MRI(3D-T1WI)への適用
3.2. 拡散MRIへの適用


### 1. 目的

  • MRtrixを用いたギブズのリンギングアーチファクト(Gibbs ringing)の除去

2. コマンド

ギブズのリンギングアーチファクト(Gibbs ringing)を除去するには、MRtrixmrdegibbsを用いる。

mrdegibbsのヘルプは、以下の通り。

クリックして展開
SYNOPSIS

     Remove Gibbs Ringing Artifacts

USAGE

     mrdegibbs [ options ] in out

        in           the input image.

        out          the output image.


DESCRIPTION

     This application attempts to remove Gibbs ringing artefacts from MRI
     images using the method of local subvoxel-shifts proposed by Kellner et
     al. (see reference below for details).

     This command is designed to run on data directly after it has been
     reconstructed by the scanner, before any interpolation of any kind has
     taken place. You should not run this command after any form of motion
     correction (e.g. not after dwifslpreproc). Similarly, if you intend
     running dwidenoise, you should run denoising before this command to not
     alter the noise structure, which would impact on dwidenoise's performance.

     Note that this method is designed to work on images acquired with full
     k-space coverage. Running this method on partial Fourier ('half-scan')
     data may lead to suboptimal and/or biased results, as noted in the
     original reference below. There is currently no means of dealing with
     this; users should exercise caution when using this method on partial
     Fourier data, and inspect its output for any obvious artefacts. 

OPTIONS

  -axes list
     select the slice axes (default: 0,1 - i.e. x-y).

  -nshifts value
     discretization of subpixel spacing (default: 20).

  -minW value
     left border of window used for TV computation (default: 1).

  -maxW value
     right border of window used for TV computation (default: 3).

Data type options

  -datatype spec
     specify output image data type. Valid choices are: float32, float32le,
     float32be, float64, float64le, float64be, int64, uint64, int64le,
     uint64le, int64be, uint64be, int32, uint32, int32le, uint32le, int32be,
     uint32be, int16, uint16, int16le, uint16le, int16be, uint16be, cfloat32,
     cfloat32le, cfloat32be, cfloat64, cfloat64le, cfloat64be, int8, uint8,
     bit.

Standard options

  -info
     display information messages.

  -quiet
     do not display information messages or progress status; alternatively,
     this can be achieved by setting the MRTRIX_QUIET environment variable to a
     non-empty string.

  -debug
     display debugging messages.

  -force
     force overwrite of output files (caution: using the same file as input and
     output might cause unexpected behaviour).

  -nthreads number
     use this number of threads in multi-threaded applications (set to 0 to
     disable multi-threading).

  -config key value  (multiple uses permitted)
     temporarily set the value of an MRtrix config file entry.

  -help
     display this information page and exit.

  -version
     display version information and exit.

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

mrdegibbs <入力画像> <出力画像> -axes 0,1  # Axial収集
mrdegibbs <入力画像> <出力画像> -axes 0,2  # Coronal収集
mrdegibbs <入力画像> <出力画像> -axes 1,2  # Sagittal収集

3. 使用例

3.1. 構造MRI(3D-T1WI)への適用

3D-T1WI(T1w.nii.gz)を入力としてmrdegibbsを実行する。

mrdegibbs T1w.nii.gz T1w_unringed.nii.gz -axes 0,1

実行した結果は、以下の通り。

3.2. 拡散MRIへの適用

mrdegibbsは、拡散MRIにも適用できる。ただし以下の条件がある。

  • dwidenoiseでノイズを除去した後に実行
  • dwifslpreprocのような歪み・頭の動き補正をする前に実行

dwidenoiseでノイズ除去された拡散強調像(DWI_denoised.nii.gz)を入力としてmrdegibbsを実行するには、以下のコマンドを実行する。

mrdegibbs DWI_denoised.nii.gz DWI_denoised_unringed.nii.gz -axes 0,1

実行した結果は、以下の通り。

【FSL】大脳基底核のセグメンテーション


1. 目的
2. コマンド
3. 使用例
3.1. 頭蓋除去をしていない場合
3.2. 頭蓋除去を既にしている場合
4. 結果
5. おまけ


1. 目的

  • 大脳基底核のセグメンテーション

2. コマンド

FSLrun_first_allを用いる。

run_first_allのヘルプは、以下。

Usage: run_first_all [options] -i <input_image> -o <output_image>

Optional arguments:
  -m <method>      : method must be one of auto, fast, none or a (numerical) threshold value
  -b               : input is already brain extracted
  -s <name>        : run only on one specified structure (e.g. L_Hipp) or a comma separated list (no spaces)
  -a <img2std.mat> : use affine matrix (do not re-run registration)
  -3               : use 3-stage affine registration (only currently for hippocampus)
  -d               : do not cleanup image output files (useful for debugging)
  -v               : verbose output
  -h               : display this help message

e.g.:  run_first_all -i im1 -o output_name 

基本的な使い方は、次の通り。

頭蓋除去をしていない3D-T1WIを入力する場合。

run_first_all -i <入力画像(3D-T1WI)> -o <出力画像(の接頭辞)>

頭蓋除去済みの3D-T1WIを入力する場合。

run_first_all -i <入力画像(3D-T1WI)> -o <出力画像(の接頭辞)> -b

3. 使用例

3.1. 頭蓋除去をしていない場合

頭蓋除去をしていない3D-T1WI(T1w.nii.gz)に対して、run_first_allをかけるには、次のようにコマンドを実行する。

ここでは、出力画像の接頭辞(オプション:-o)を「output」とした。

run_first_all -i T1w.nii.gz -o output

3.2. 頭蓋除去を既にしている場合

まず、頭蓋除去済みの3D-T1WI(T1_skull_stripped.nii.gz)を用意する。

頭蓋除去のやり方は、以下の記事を参考にするとよい。

頭蓋除去した画像(T1_skull_stripped.nii.gz)に対して、run_first_allコマンドを実行する。

ここでは、出力画像の接頭辞(オプション:-o)を「output」とし、頭蓋除去済みの脳を入力していることを示すオプション-bを付けている。

run_first_all -i T1_skull_stripped.nii.gz -o output -b

4. 結果

処理が完了すると、次のファイルが出力される。

  • output_name_all_fast_firstseg.nii.gz:大脳基底核がセグメンテーションされた画像(3D画像)
  • output_name_all_fast_origsegs.nii.gz:大脳基底核の各セグメンテーション画像(4D画像)
  • output_name_first.vtk:セグメンテーションをする際のメッシュ。FSLViewの3Dモードで見ることができる。
  • output_name_first.bvars:パラメータファイル。

大脳基底核のセグメント(output_all_fast_firstseg.nii.gz)と頭蓋除去した画像(T1_skull_stripped.nii.gz)を、重ね合わせてみる。

fsleyes T1_skull_stripped.nii.gz output_all_fast_firstseg.nii.gz -cm random

run_first_allでは、大脳基底核を次の領域にセグメント(区域分け)する。

Label Index 領域
10 Left-Thalamus-Proper
11 Left-Caudate
12 Left-Putamen
13 Left-Pallidum
16 Brain-Stem/ 4th Ventricle
17 Left-Hippocampus
18 Left-Amygdala
26 Left-Accumbens-area
49 Right-Thalamus-Proper
50 Right-Caudate
51 Right-Putamen
52 Right-Pallidum
53 Right-Hippocampus
54 Right-Amygdala
58 Right-Accumbens-area

5. おまけ

複数の被験者のセグメンテーション結果をQCしたい場合、first_roi_slicesdirコマンドを用いると便利である。

基本な使い方は、以下。

first_roi_slicesdir <入力画像(3D-T1WI)のリスト> <セグメンテーション画像のリスト>

例えば、次のような被験者Subj001, Subj002, Subj003がいた場合。

.
├── Subj001_T1_skull_stripped.nii.gz
├── Subj001_output_all_fast_firstseg.nii.gz
├── Subj002_T1_skull_stripped.nii.gz
├── Subj002_output_all_fast_firstseg.nii.gz
├── Subj003_T1_skull_stripped.nii.gz
└── Subj003_output_all_fast_firstseg.nii.gz

first_roi_slicesdirを、次のように実行する。

first_roi_slicesdir *_t1.nii.gz *_all_fast_firstseg.nii.gz

処理が完了すると、「slicesdir」フォルダが生成される。

slicesdir/
├── Subj001_T1_skull_stripped_t1grot1_to_Subj001_output_all_fast_firstseglbgrot1.png
├── Subj002_T1_skull_stripped_t1grot2_to_Subj002_output_all_fast_firstseglbgrot2.png
├── Subj003_T1_skull_stripped_t1grot3_to_Subj003_output_all_fast_firstseglbgrot3.png
├── grota.png
├── grotb.png
├── grotc.png
├── grotd.png
├── grote.png
├── grotf.png
├── grotg.png
├── groth.png
├── groti.png
└── index.html

slicesdirフォルダの「index.html」を開くと、結果が見れる。

【FreeSurfer】FreeSurferを用いた頭蓋除去 ~Skull-stripping~


1. 目的
2. コマンド
3. 使用例
4. 結果


1. 目的

脳画像解析で、対象となる領域は脳実質でありその他の骨・筋・脂肪・眼球等の組織は、解析する上でノイズとなる。そのため、解析精度を高めるには頭蓋を除去することが重要である。

この記事では、FreeSurferを用いて頭蓋を除去する方法を解説する。

  • FreeSurferを用いた頭蓋除去

2. コマンド

ここでは、FreeSurferの関数であるrecon-allコマンドを用いて、頭蓋を除去する。

基本的な使い方は、次の通り。

# FreeSurfer作業ディレクトリを指定
export SUBJECTS_DIR=<PATH>

# Autorecon1を実行
recon-all -subjid <被験者ID (任意)> -i <入力画像 (3D-T1WI)> -autorecon1

Autorecon1では、以下の処理を実行している。

  1. 動き補正(同一被験者の3D-T1WIが二つある場合には、平均画像を生成)
  2. Talairach変換
  3. 信号ムラ(バイアス)の補正
  4. 信号値の正規化
  5. 頭蓋除去

3. 使用例

頭蓋除去前の3D-T1WI(T1w.nii.gz)に対して、recon-all -autorecon1コマンドを実行する。ここでは、被験者ID(-subjid)を「Subj001」にした。これは、自由に変えてもよい。

# FreeSurfer作業ディレクトリを指定
export SUBJECTS_DIR=$PWD

# Autorecon1を実行
recon-all -subjid Subj001 -i T1w.nii.gz -autorecon1

# Autorecon1で作成した頭蓋除去済みの脳マスク画像(brainmask)をNIfTIに変換
# ただし、元画像(T1w.nii.gz)と同じヘッダー情報を参考に(--likeオプション)
mri_convert Subj001/mri/brainmask.mgz brainmask.nii.gz --like T1w.nii.gz

# 元のT1WIにbrainmaskを適用
fslmaths T1w.nii.gz -mas brainmask.nii.gz T1_skull_stripped.nii.gz

recon-allでの処理が完了すると、被験者IDディレクトリ(Subj001)が生成され、このディレクトリ内に処理された結果が保存される。

Subj001/
├── label
├── mri
│   ├── T1.mgz
│   ├── brainmask.auto.mgz
│   ├── brainmask.mgz
│   ├── mri_nu_correct.mni.log
│   ├── mri_nu_correct.mni.log.bak
│   ├── nu.mgz
│   ├── orig
│   │   └── 001.mgz
│   ├── orig.mgz
│   ├── orig_nu.mgz
│   ├── rawavg.mgz
│   ├── talairach_with_skull.log
│   └── transforms
│       ├── bak
│       ├── talairach.auto.xfm
│       ├── talairach.auto.xfm.lta
│       ├── talairach.xfm
│       ├── talairach.xfm.lta
│       ├── talairach_avi.log
│       ├── talairach_avi_QA.log
│       ├── talairach_with_skull.lta
│       └── talsrcimg_to_711-2C_as_mni_average_305_t4_vox2vox.txt
├── scripts
│   ├── build-stamp.txt
│   ├── lastcall.build-stamp.txt
│   ├── patchdir.txt
│   ├── recon-all-status.log
│   ├── recon-all.cmd
│   ├── recon-all.done
│   ├── recon-all.env
│   ├── recon-all.local-copy
│   ├── recon-all.log
│   ├── recon-config.yaml
│   └── unknown-args.txt
├── stats
├── surf
├── tmp
├── touch
│   ├── conform.touch
│   ├── inorm1.touch
│   ├── nu.touch
│   ├── skull.lta.touch
│   ├── skull_strip.touch
│   └── talairach.touch
└── trash

頭蓋除去された画像は、「Subj001/mri/brainmask.mgz」にある。brainmask.mgz画像は、頭蓋除去済みの脳画像であるが、信号値の階調が0~225に圧縮(正規化)されてしまう。そこで、このbrainmask.mgzをマスク画像として扱い、元の3D-T1WIのマスク処理(マスキング)として適用する。

結果

処理前後の画像は次の通り。

【FSL】FSLを用いた頭蓋除去 ~Skull-stripping~


1. 目的
2. コマンド
3. 使用例
4. 結果


1. 目的

脳画像解析で、対象となる領域は脳実質でありその他の骨・筋・脂肪・眼球等の組織は、解析する上でノイズとなる。そのため、解析精度を高めるには頭蓋を除去することが重要である。

この記事では、FSLを用いて頭蓋を除去する方法を解説する。

  • FSLを用いた頭蓋除去

2. コマンド

ここでは、FSLの関数であるBETコマンドを用いて、頭蓋を除去する。

betコマンドのヘルプは次の通り。


Usage:    bet <input> <output> [options]

Main bet2 options:
  -o          generate brain surface outline overlaid onto original image
  -m          generate binary brain mask
  -s          generate approximate skull image
  -n          don't generate segmented brain image output
  -f <f>      fractional intensity threshold (0->1); default=0.5; smaller values give larger brain outline estimates
  -g <g>      vertical gradient in fractional intensity threshold (-1->1); default=0; positive values give larger brain outline at bottom, smaller at top
  -r <r>      head radius (mm not voxels); initial surface sphere is set to half of this
  -c <x y z>  centre-of-gravity (voxels not mm) of initial mesh surface.
  -t          apply thresholding to segmented brain image and mask
  -e          generates brain surface as mesh in .vtk format

Variations on default bet2 functionality (mutually exclusive options):
  (default)   just run bet2
  -R          robust brain centre estimation (iterates BET several times)
  -S          eye & optic nerve cleanup (can be useful in SIENA - disables -o option)
  -B          bias field & neck cleanup (can be useful in SIENA)
  -Z          improve BET if FOV is very small in Z (by temporarily padding end slices)
  -F          apply to 4D FMRI data (uses -f 0.3 and dilates brain mask slightly)
  -A          run bet2 and then betsurf to get additional skull and scalp surfaces (includes registrations)
  -A2 <T2>    as with -A, when also feeding in non-brain-extracted T2 (includes registrations)

Miscellaneous options:
  -v          verbose (switch on diagnostic messages)
  -h          display this help, then exits
  -d          debug (don't delete temporary intermediate images)

主要なオプションの役割は、次の通り。

  • -o:脳表の輪郭画像を生成
  • -m:バイナリ―のマスク画像を生成
  • -s:骨画像を生成
  • -n:セグメントした画像を生成しない
  • -f :信号値に対するしきい値(0-1, 初期値:0.5). 小さいほど脳が大きく出力される.
  • -g :信号値しきい値における垂直方向の勾配. 大きくなるほど脳幹付近の脳が大きくなり、頭頂部の脳が小さくなる.
  • -t:セグメントされた脳とマスク画像にしきい値処理を適用
  • -R:堅牢性の高い脳の中心推定(BETを何回か繰り返す)
  • -S :眼球と視神経を除去
  • -B :バイアスフィールドの補正と首の除去
  • -Z:z方向(スライス方向)のFOVがかなり小さい場合にBETの精度を上げる(一時的にスライスをパディング)
    • -F:4DのfMRIデータに適用(”-f 0.3″として少し脳マスクを大きくする)

基本的な使い方は、次の通り。

bet <入力画像> <出力画像> [オプション]

オプションについては、出力画像の結果を見ながら調節するとよい。

3. 使用例

頭蓋除去前の3D-T1WI(T1w.nii.gz)に対して、betコマンドを実行する。

bet T1w.nii.gz T1_skull_stripped.nii.gz -f 0.3 -R -S -B

処理が完了すると、頭蓋除去された画像(T1_skull_stripped.nii.gz)とそのマスク画像(T1_skull_stripped_mask.nii.gz)が生成される。

ls  # カレントディレクトリのファイルを確認

Out:

    T1_skull_stripped.nii.gz  T1_skull_stripped_mask.nii.gz  T1w.nii.gz

4. 結果

処理前後の画像は次の通り。

【FreeSurfer】FreeSurferを用いたRF/B1バイアス(信号ムラ)補正


1. 目的
2. コマンド
3. 使用例
4. 結果


1. 目的

  • RF/B1バイアス(信号ムラ)補正

2. コマンド

ここでは、FreeSurferの関数であるAntsN4BiasFieldCorrectionFsコマンドを用いて、信号ムラ補正をする。

AntsN4BiasFieldCorrectionFsは、ANTsの関数である”AntsN4BiasFieldCorrection”を基にしたものであり、補正アルゴリズムとしては、ANTs由来のものである。

AntsN4BiasFieldCorrectionFsコマンドのヘルプは次の通り。

				Help

NAME
	AntsN4BiasFieldCorrectionFs

SYNOPSIS
	AntsN4BiasFieldCorrectionFs [options] -i <invol> -o <outvol>

DESCRIPTION
	Runs N4 (nonparameteric, nonuniform normalization) retrospective bias 
	correction on an image. This programs wraps the 
	AntsN4BiasFieldCorrection utility available in the ANTs package (see 
	http://stnava.github.io/ANTs).

REQUIRED FLAGGED ARGUMENTS
	-i, --input invol
		input volume file

	-o, --output outvol
		corrected volume file

OPTIONAL FLAGGED ARGUMENTS
	--shrink
		resample factor to decrease computation time (default is 4)

基本的な使い方は、次の通り。

AntsN4BiasFieldCorrectionFs -i <入力画像> -o <出力画像>

3. 使用例

補正前の3D-T1WI(T1w.nii.gz)に対して、信号ムラ補正をする。

AntsN4BiasFieldCorrectionFs -i T1w.nii.gz -o T1w_biascorrected.nii.gz

処理が完了すると、補正後の3D-T1WI(T1w_biascorrected.nii.gz)が出力される。

ls  # カレントディレクトリのファイルを確認
    T1w.nii.gz    T1w_biascorrected.nii.gz

4. 結果

補正前(上)と補正後(下)を比較すると次の通り。

補正前では頭頂部で比較的低信号、深部灰白質で比較的高信号だったのが、補正後に均一になっている。

【FSL】 FSLを用いたRF/B1バイアス(信号ムラ)補正とセグメンテーション


1. 目的
2. コマンド
3. 使用例
4. 結果
4.1. バイアス(信号ムラ)補正
4.2. ボクセルあたりの部分容積(存在割合)
4.3. 脳組織のセグメンテーション
4.4. 確率マップ


1. 目的

  • RF/B1バイアス(信号ムラ)補正
  • 脳組織のセグメンテーション(Segmentation)

2. コマンド

FSLfastコマンドを用いる。

fastコマンドヘルプは次の通り。

Usage: 
fast [options] file(s)

Optional arguments (You may optionally specify one or more of):
	-n,--class	number of tissue-type classes; default=3
	-I,--iter	number of main-loop iterations during bias-field removal; default=4
	-l,--lowpass	bias field smoothing extent (FWHM) in mm; default=20
	-t,--type	type of image 1=T1, 2=T2, 3=PD; default=T1
	-f,--fHard	initial segmentation spatial smoothness (during bias field estimation); default=0.02
	-g,--segments	outputs a separate binary image for each tissue type
	-a <standard2input.mat> initialise using priors; you must supply a FLIRT transform
	-A <prior1> <prior2> <prior3>    alternative prior images
	--nopve	turn off PVE (partial volume estimation)
	-b		output estimated bias field
	-B		output bias-corrected image
	-N,--nobias	do not remove bias field
	-S,--channels	number of input images (channels); default 1
	-o,--out	output basename
	-P,--Prior	use priors throughout; you must also set the -a option
	-W,--init	number of segmentation-initialisation iterations; default=15
	-R,--mixel	spatial smoothness for mixeltype; default=0.3
	-O,--fixed	number of main-loop iterations after bias-field removal; default=4
	-H,--Hyper	segmentation spatial smoothness; default=0.1
	-v,--verbose	switch on diagnostic messages
	-h,--help	display this message
	-s,--manualseg <filename> Filename containing intensities
	-p		outputs individual probability maps

基本的な使い方は、次の通り。

入力画像として、T1WI, T2WI, PDWIを用いることができる。

fast <入力画像> T1w.nii.gz

3. 使用例

まず前処理として、処理前の3D-T1WI(T1w.nii.gz)に対して、頭蓋除去(skull-stripping)をする。頭蓋除去のやり方は、以下を参考に。

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

頭蓋除去した画像(T1_skull_stripped.nii.gz)に対して、fastコマンドを実行する。

fast -t 1 -g -B -b -p -o output T1_skull_stripped.nii.gz

ここで、使用したオプションは次の通り。

  • -t:入力画像の種類。ここでは、T1WIであるため「-t 1」と設定。
  • -g:セグメント後の画像をバイナリー画像として出力
  • -b:バイアスフィールドを出力
  • -B:バイアス補正後の画像を出力
  • -p:セグメントした各確率マップを出力
  • -o:出力画像の接頭辞を指定。ここでは「output」とした。

4. 結果

処理が完了すると、次のような画像が出力される。

4.1. バイアス(信号ムラ)補正

  • output_bias.nii.gz:バイアスフィールド画像
  • output_restore.nii.gz:バイアス補正後の画像

4.2. ボクセルあたりの部分容積(存在割合)

  • output_pve_0.nii.gz:脳脊髄液の部分容積画像
  • output_pve_1.nii.gz:灰白質の部分容積画像
  • output_pve_2.nii.gz:白質の部分容積画像

4.3. 脳組織のセグメンテーション

  • output_seg_0.nii.gz:脳脊髄液の2値(バイナリー)画像
  • output_seg_1.nii.gz:灰白質の2値画像
  • output_seg_2.nii.gz:白質の2値画像
  • output_seg.nii.gz:脳脊髄液(1)・灰白質(2)・白質(3)がセグメントされた画像

4.4. 確率マップ

  • output_prob_0.nii.gz:脳脊髄液の確率マップ
  • output_prob_1.nii.gz:灰白質の確率マップ
  • output_prob_2.nii.gz:白質の確率マップ

【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

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

Anacondaに頼らない、pipとvenvを用いたPython環境の構築

最近、Pythonに触れることが多くなってきました。
その中で、環境構築についていろいろ学んできました。

Pythonの参考書の多くは”Anacondaで環境構築しましょう”と書いてあります。
しかし、Anacondaはセットアップファイルだけで4GBもあります。
また、自分のシステムに既に入っているPythonとの相互関係も最初の頃はよくわからなくなります。

Anacondaを横においておくと、Pythonには、パッケージマネージャーとして、”pip” というものがあります。
これも若干クセがあるので、いくつかおさえておくべきことがあります。

さらに、Pythonは”venv”というパッケージを使うことで、仮想環境を簡単に構築できます。
このvenvについて把握すると、Anacondaなどのことも理解しやすくなります。

ということで、私なりに理解したことをここでまとめていきたいと思います。
なお、ここではすべてPython3環境を意識していきます。pipはmacOSやUbuntuでは全部Python3になっています。DebianではPython2のようですが、最近、Debianを使っていないのでよくわかりません。(man pip に書いてある情報から記載しただけです)

現時点での私のおすすめは、
「基本、–userをつけてpipでインストール。試験的に試したかったらvenvで仮想環境内で構築」です。

概要は以下になります。

  1. pip
    1.1 どのpipを使っているかの確認
    1.2 システムへのインストールと個別ユーザーへのインストール
  2. venv
    2.1 仮想環境の構築
    2.2 仮想環境の有効化
    2.3 仮想環境の無効化

続きを読む

macOSでのSPM12のコンパイル方法

ある方から、Apple M1のmacでSPMを起動しようとするとspm_check_installation(‘basic’)でエラーが出て起動しないという相談を受けました。

コンパイルしたら問題は解決しました。コンパイル方法を共有します。

ただし、その後、SPMのMLでこのディスカッションに乗ってみたところ、コンパイルは不要だよということも教えていただきました。なので、コンパイルに挑戦してみたい人向けと思ってください。(普通は不要です)

続きを読む

格安パルスオキシメーターは使えるのか?

久しぶりに脳画像以外のネタを。

COVID-19が猛威をふるう中、動脈血酸素飽和度 (SpO2) を測定できるパルスオキシメーターの需要が増えています。
万が一自宅療養になる時などにそなえてパルスオキシメーターが家にあるといいなと思いましたが、
一般的に購入できる1万円未満のパルスオキシメーターのレビューを見ると、みなさん、様々なことを書いていて、判断に困るなぁと思いました。

そこで、実際どうなのかと思い、自ら人柱になって、自腹で購入して実験してみました。

続きを読む