【FSL/MRtrix】FSL/MRtrixを用いたしきい値処理とマスク画像の作成


1. 目的
2. FSLを用いる場合
2.1. コマンド
2.2. ノイズ除去(デノイズ)
2.3. 複数のラベルから1部のラベルを抽出
3. MRtrixを用いる場合
3.1. コマンド
3.2. ノイズ除去(デノイズ)
3.3. 複数のラベルから1部のラベルを抽出


1. 目的

  • FSL/MRtrixを用いたしきい値処理とマスク画像の作成

2. FSLを用いる場合

2.1. コマンド

FSLfslmathsを用いる。fslmathsは、画像の四則演算からしきい値処理、フィルタリングなど基本的な画像処理を実行することができるコマンドである。

fslmathsのヘルプは、次の通り。

クリックして展開
Usage: fslmaths [-dt <datatype>] <first_input> [operations and inputs] <output> [-odt <datatype>]

Datatype information:
 -dt sets the datatype used internally for calculations (default float for all except double images)
 -odt sets the output datatype ( default is float )
 Possible datatypes are: char short int float double input
 "input" will set the datatype to that of the original image

Binary operations:
  (some inputs can be either an image or a number)
 -add   : add following input to current image
 -sub   : subtract following input from current image
 -mul   : multiply current image by following input
 -div   : divide current image by following input
 -rem   : modulus remainder - divide current image by following input and take remainder
 -mas   : use (following image>0) to mask current image
 -thr   : use following number to threshold current image (zero anything below the number)
 -thrp  : use following percentage (0-100) of ROBUST RANGE to threshold current image (zero anything below the number)
 -thrP  : use following percentage (0-100) of ROBUST RANGE of non-zero voxels and threshold below
 -uthr  : use following number to upper-threshold current image (zero anything above the number)
 -uthrp : use following percentage (0-100) of ROBUST RANGE to upper-threshold current image (zero anything above the number)
 -uthrP : use following percentage (0-100) of ROBUST RANGE of non-zero voxels and threshold above
 -max   : take maximum of following input and current image
 -min   : take minimum of following input and current image
 -seed  : seed random number generator with following number
 -restart : replace the current image with input for future processing operations
 -save : save the current working image to the input filename

Basic unary operations:
 -exp   : exponential
 -log   : natural logarithm
 -sin   : sine function
 -cos   : cosine function
 -tan   : tangent function
 -asin  : arc sine function
 -acos  : arc cosine function
 -atan  : arc tangent function
 -sqr   : square
 -sqrt  : square root
 -recip : reciprocal (1/current image)
 -abs   : absolute value
 -bin   : use (current image>0) to binarise
 -binv  : binarise and invert (binarisation and logical inversion)
 -fillh : fill holes in a binary mask (holes are internal - i.e. do not touch the edge of the FOV)
 -fillh26 : fill holes using 26 connectivity
 -index : replace each nonzero voxel with a unique (subject to wrapping) index number
 -grid <value> <spacing> : add a 3D grid of intensity <value> with grid spacing <spacing>
 -edge  : edge strength
 -tfce <H> <E> <connectivity>: enhance with TFCE, e.g. -tfce 2 0.5 6 (maybe change 6 to 26 for skeletons)
 -tfceS <H> <E> <connectivity> <X> <Y> <Z> <tfce_thresh>: show support area for voxel (X,Y,Z)
 -nan   : replace NaNs (improper numbers) with 0
 -nanm  : make NaN (improper number) mask with 1 for NaN voxels, 0 otherwise
 -rand  : add uniform noise (range 0:1)
 -randn : add Gaussian noise (mean=0 sigma=1)
 -inm <mean> :  (-i i ip.c) intensity normalisation (per 3D volume mean)
 -ing <mean> :  (-I i ip.c) intensity normalisation, global 4D mean)
 -range : set the output calmin/max to full data range

Matrix operations:
 -tensor_decomp : convert a 4D (6-timepoint )tensor image into L1,2,3,FA,MD,MO,V1,2,3 (remaining image in pipeline is FA)

Kernel operations (set BEFORE filtering operation if desired):
 -kernel 3D : 3x3x3 box centered on target voxel (set as default kernel)
 -kernel 2D : 3x3x1 box centered on target voxel
 -kernel box    <size>     : all voxels in a cube of width <size> mm centered on target voxel
 -kernel boxv   <size>     : all voxels in a cube of width <size> voxels centered on target voxel, CAUTION: size should be an odd number
 -kernel boxv3  <X> <Y> <Z>: all voxels in a cuboid of dimensions X x Y x Z centered on target voxel, CAUTION: size should be an odd number
 -kernel gauss  <sigma>    : gaussian kernel (sigma in mm, not voxels)
 -kernel sphere <size>     : all voxels in a sphere of radius <size> mm centered on target voxel
 -kernel file   <filename> : use external file as kernel

Spatial Filtering operations: N.B. all options apart from -s use the default kernel or that _previously_ specified by -kernel
 -dilM    : Mean Dilation of non-zero voxels
 -dilD    : Modal Dilation of non-zero voxels
 -dilF    : Maximum filtering of all voxels
 -dilall  : Apply -dilM repeatedly until the entire FOV is covered
 -ero     : Erode by zeroing non-zero voxels when zero voxels found in kernel
 -eroF    : Minimum filtering of all voxels
 -fmedian : Median Filtering 
 -fmean   : Mean filtering, kernel weighted (conventionally used with gauss kernel)
 -fmeanu  : Mean filtering, kernel weighted, un-normalised (gives edge effects)
 -s <sigma> : create a gauss kernel of sigma mm and perform mean filtering
 -subsamp2  : downsamples image by a factor of 2 (keeping new voxels centred on old)
 -subsamp2offc  : downsamples image by a factor of 2 (non-centred)

Dimensionality reduction operations:
  (the "T" can be replaced by X, Y or Z to collapse across a different dimension)
 -Tmean   : mean across time
 -Tstd    : standard deviation across time
 -Tmax    : max across time
 -Tmaxn   : time index of max across time
 -Tmin    : min across time
 -Tmedian : median across time
 -Tperc <percentage> : nth percentile (0-100) of FULL RANGE across time
 -Tar1    : temporal AR(1) coefficient (use -odt float and probably demean first)

Basic statistical operations:
 -pval    : Nonparametric uncorrected P-value, assuming timepoints are the permutations; first timepoint is actual (unpermuted) stats image
 -pval0   : Same as -pval, but treat zeros as missing data
 -cpval   : Same as -pval, but gives FWE corrected P-values
 -ztop    : Convert Z-stat to (uncorrected) P
 -ptoz    : Convert (uncorrected) P to Z
 -rank    : Convert data to ranks (over T dim)
 -ranknorm: Transform to Normal dist via ranks

Multi-argument operations:
 -roi <xmin> <xsize> <ymin> <ysize> <zmin> <zsize> <tmin> <tsize> : zero outside roi (using voxel coordinates). Inputting -1 for a size will set it to the full image extent for that dimension.
 -bptf  <hp_sigma> <lp_sigma> : (-t in ip.c) Bandpass temporal filtering; nonlinear highpass and Gaussian linear lowpass (with sigmas in volumes, not seconds); set either sigma<0 to skip that filter
 -roc <AROC-thresh> <outfile> [4Dnoiseonly] <truth> : take (normally binary) truth and test current image in ROC analysis against truth. <AROC-thresh> is usually 0.05 and is limit of Area-under-ROC measure FP axis. <outfile> is a text file of the ROC curve (triplets of values: FP TP threshold). If the truth image contains negative voxels these get excluded from all calculations. If <AROC-thresh> is positive then the [4Dnoiseonly] option needs to be set, and the FP rate is determined from this noise-only data, and is set to be the fraction of timepoints where any FP (anywhere) is seen, as found in the noise-only 4d-dataset. This is then controlling the FWE rate. If <AROC-thresh> is negative the FP rate is calculated from the zero-value parts of the <truth> image, this time averaging voxelwise FP rate over all timepoints. In both cases the TP rate is the average fraction of truth=positive voxels correctly found.

Combining 4D and 3D images:
 If you apply a Binary operation (one that takes the current image and a new image together), when one is 3D and the other is 4D,
 the 3D image is cloned temporally to match the temporal dimensions of the 4D image.

e.g. fslmaths inputVolume -add inputVolume2 output_volume
     fslmaths inputVolume -add 2.5 output_volume
     fslmaths inputVolume -add 2.5 -mul inputVolume2 output_volume

     fslmaths 4D_inputVolume -Tmean -mul -1 -add 4D_inputVolume demeaned_4D_inputVolume

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

fslmaths  <入力画像1> [演算子あるいは入力画像] <出力画像> 

2.2. ノイズ除去(デノイズ)

ここでは、特にしきい値処理で用いる-thr-uthrオプション、さらにバイナリーマスク作成に必要な-binオプションを例にfslmathsコマンドの使い方を解説する。

例えば、拡散MRI(b=0, SE-EPI, DWI_b0.nii.gz)に対して、二値化処理し脳マスク画像を生成する場合、以下のようなコマンドになる。

fslmaths DWI_b0.nii.gz -bin DWI_b0_mask.nii.gz

生成した脳マスク画像(緑)と拡散MRI(b=0, SE-EPI, DWI_b0.nii.gz)を重ね合わせてみる。脳以外の領域に至るまでマスキングしていることが分かる。

脳周囲のノイズ信号値を確認すると、0~30程度であった。

そこで、信号値30以下をカットするようにしきい値処理をするために、-thrオプションを用いる。

fslmaths DWI_b0.nii.gz -thr 30 -bin DWI_b0_mask_thr30.nii.gz

しきい値処理をして生成した脳マスク画像(緑)と拡散MRI(b=0, SE-EPI, DWI_b0.nii.gz)を重ね合わせてみる。ノイズ部分のマスキングが解消されていることが分かる。

2.3. 複数のラベルから1部のラベルを抽出

以下のような、CSF/GM/WMのラベル(CSF: 1, GM: 2, WM: 3)があったとする。

この内、GMのみを抽出したい場合、下限値-thrおよび上限値-uthr共に信号値2になるように設定すればよい。

fslmaths CSF_GM_WM_seg.nii.gz -thr 2 -uthr 2 GM.nii.gz

CSF/GM/WMのラベルからGMのラベルのみが抽出される。

3. MRtrixを用いる場合

3.1. コマンド

MRtrixmrthresholdを用いる。mrthresholdは、画像のしきい値処理に用いるコマンドである。

mrthresholdのヘルプは、次の通り。

クリックして展開
USAGE

     mrthreshold [ options ] input [ output ]

        input        the input image to be thresholded

        output       the (optional) output binary image mask


DESCRIPTION

     The threshold value to be applied can be determined in one of a number of
     ways:

     - If no relevant command-line option is used, the command will
     automatically determine an optimal threshold;

     - The -abs option provides the threshold value explicitly;

     - The -percentile, -top and -bottom options enable more fine-grained
     control over how the threshold value is determined.

     The -mask option only influences those image values that contribute toward
     the determination of the threshold value; once the threshold is
     determined, it is applied to the entire image, irrespective of use of the
     -mask option. If you wish for the voxels outside of the specified mask to
     additionally be excluded from the output mask, this can be achieved by
     providing the -out_masked option.

     The four operators available through the "-comparison" option ("lt", "le",
     "ge" and "gt") correspond to "less-than" (<), "less-than-or-equal" (<=),
     "greater-than-or-equal" (>=) and "greater-than" (>). This offers
     fine-grained control over how the thresholding operation will behave in
     the presence of values equivalent to the threshold. By default, the
     command will select voxels with values greater than or equal to the
     determined threshold ("ge"); unless the -bottom option is used, in which
     case after a threshold is determined from the relevant lowest-valued image
     voxels, those voxels with values less than or equal to that threshold
     ("le") are selected. This provides more fine-grained control than the
     -invert option; the latter is provided for backwards compatibility, but is
     equivalent to selection of the opposite comparison within this selection.

     If no output image path is specified, the command will instead write to
     standard output the determined threshold value.

Threshold determination mechanisms

  -abs value
     specify threshold value as absolute intensity

  -percentile value
     determine threshold based on some percentile of the image intensity
     distribution

  -top count
     determine threshold that will result in selection of some number of
     top-valued voxels

  -bottom count
     determine & apply threshold resulting in selection of some number of
     bottom-valued voxels (note: implies threshold application operator of "le"
     unless otherwise specified)

Threshold determination modifiers

  -allvolumes
     compute a single threshold for all image volumes, rather than an
     individual threshold per volume

  -ignorezero
     ignore zero-valued input values during threshold determination

  -mask image
     compute the threshold based only on values within an input mask image

Threshold application modifiers

  -comparison choice
     comparison operator to use when applying the threshold; options are:
     lt,le,ge,gt (default = "le" for -bottom; "ge" otherwise)

  -invert
     invert the output binary mask (equivalent to flipping the operator;
     provided for backwards compatibility)

  -out_masked
     mask the output image based on the provided input mask image

  -nan
     set voxels that fail the threshold to NaN rather than zero (output image
     will be floating-point rather than binary)

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.

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

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

3.2. ノイズ除去(デノイズ)

例えば、拡散MRI(b=0, SE-EPI, DWI_b0.nii.gz)に対して、二値化処理し脳マスク画像を生成する場合、以下のようなコマンドになる。

ここで、-absはしきい値を設定するオプションであり、-comparisonはしきい値に対してどのような操作を実行するのかを指定するオプションである。例えば、-comparisonでは、の4種類(“lt”, “le”, “ge”, “gt”)の操作ができ、それぞれ“less-than” (<), “less-than-or-equal” (<=), “greater-than-or-equal” (>=), “greater-than” (>)を意味する。

mrthreshold -abs 0 -comparison gt DWI_b0.nii.gz DWI_b0_mask.nii.gz

生成した脳マスク画像(緑)と拡散MRI(b=0, SE-EPI, DWI_b0.nii.gz)を重ね合わせてみる。脳以外の領域に至るまでマスキングしていることが分かる。

脳周囲のノイズ信号値を確認すると、0~30程度であった。

そこで、信号値30以下をカットするようにしきい値処理をするために、-abs 30とする。

mrthreshold -abs 30 -comparison gt DWI_b0.nii.gz DWI_b0_mask_thr30.nii.gz

しきい値処理をして生成した脳マスク画像(緑)と拡散MRI(b=0, SE-EPI, DWI_b0.nii.gz)を重ね合わせてみる。ノイズ部分のマスキングが解消されていることが分かる。

3.3. 複数のラベルから1部のラベルを抽出

以下のような、CSF/GM/WMのラベル(CSF: 1, GM: 2, WM: 3)があったとする。

この内、WMのみを抽出したい場合、次のようにコマンドを実行する。

mrthreshold -abs 2 -comparison gt CSF_GM_WM_seg.nii.gz WM.nii.gz

CSF/GM/WMのラベルからWMのラベルのみが抽出される。

著者情報: 斎藤 勇哉

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

コメントを残す

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