【MRtrix】MRtrixを用いた解像度の変更 ~Upsampling~


1. 目的
2. コマンド
3. 使用例
3.1. ボクセルサイズを指定(オプション:-voxel)
3.2. スケールを指定(オプション:-scale))
3.3. ボクセルサイズを指定(オプション:-voxel))
3.4. 目的の解像度を持つ画像を指定(オプション:-template))


1. 目的

  • MRtrixを用いたアップサンプリング(Upsampling)

2. コマンド

MRtrixのmrgridを用いる。

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

クリックして展開
SYNOPSIS

     Modify the grid of an image without interpolation (cropping or padding) or
     by regridding to an image grid with modified orientation, location and or
     resolution. The image content remains in place in real world coordinates.

USAGE

     mrgrid [ options ] input operation output

        input        input image to be regridded.

        operation    the operation to be performed, one of: regrid, crop, pad.

        output       the output image.


DESCRIPTION

     - regrid: This operation performs changes of the voxel grid that require
     interpolation of the image such as changing the resolution or location and
     orientation of the voxel grid. If the image is down-sampled, the
     appropriate smoothing is automatically applied using Gaussian smoothing
     unless nearest neighbour interpolation is selected or oversample is
     changed explicitly. The resolution can only be changed for spatial
     dimensions. 

     - crop: The image extent after cropping, can be specified either manually
     for each axis dimensions, or via a mask or reference image. The image can
     be cropped to the extent of a mask. This is useful for axially-acquired
     brain images, where the image size can be reduced by a factor of 2 by
     removing the empty space on either side of the brain. Note that cropping
     does not extend the image beyond the original FOV unless explicitly
     specified (via -crop_unbound or negative -axis extent).

     - pad: Analogously to cropping, padding increases the FOV of an image
     without image interpolation. Pad and crop can be performed simultaneously
     by specifying signed specifier argument values to the -axis option.

     This command encapsulates and extends the functionality of the superseded
     commands 'mrpad', 'mrcrop' and 'mrresize'. Note the difference in -axis
     convention used for 'mrcrop' and 'mrpad' (see -axis option description).

EXAMPLE USAGES

     Crop and pad the first axis:
       $ mrgrid in.mif crop -axis 0 10,-5 out.mif
     This removes 10 voxels on the lower and pads with 5 on the upper bound,
     which is equivalent to padding with the negated specifier (mrgrid in.mif
     pad -axis 0 -10,5 out.mif).

     Right-pad the image to the number of voxels of a reference image:
       $ mrgrid in.mif pad -as ref.mif -all_axes -axis 3 0,0 out.mif -fill nan
     This pads the image on the upper bound of all axes except for the volume
     dimension. The headers of in.mif and ref.mif are ignored and the output
     image uses NAN values to fill in voxels outside the original range of
     in.mif.

     Regrid and interpolate to match the voxel grid of a reference image:
       $ mrgrid in.mif regrid -template ref.mif -scale 1,1,0.5 out.mif -fill nan
     The -template instructs to regrid in.mif to match the voxel grid of
     ref.mif (voxel size, grid orientation and voxel centres). The -scale
     option overwrites the voxel scaling factor yielding voxel sizes in the
     third dimension that are twice as coarse as those of the template image.

Regridding options (involves image interpolation, applied to spatial axes only)

  -template image
     match the input image grid (voxel spacing, image size, header
     transformation) to that of a reference image. The image resolution
     relative to the template image can be changed with one of -size, -voxel,
     -scale.

  -size dims
     define the size (number of voxels) in each spatial dimension for the
     output image. This should be specified as a comma-separated list.

  -voxel size
     define the new voxel size for the output image. This can be specified
     either as a single value to be used for all spatial dimensions, or as a
     comma-separated list of the size for each voxel dimension.

  -scale factor
     scale the image resolution by the supplied factor. This can be specified
     either as a single value to be used for all dimensions, or as a
     comma-separated list of scale factors for each dimension.

  -interp method
     set the interpolation method to use when reslicing (choices: nearest,
     linear, cubic, sinc. Default: cubic).

  -oversample factor
     set the amount of over-sampling (in the target space) to perform when
     regridding. This is particularly relevant when downsamping a
     high-resolution image to a low-resolution image, to avoid aliasing
     artefacts. This can consist of a single integer, or a comma-separated list
     of 3 integers if different oversampling factors are desired along the
     different axes. Default is determined from ratio of voxel dimensions
     (disabled for nearest-neighbour interpolation).

Pad and crop options (no image interpolation is performed, header transformation is adjusted)

  -as reference image
     pad or crop the input image on the upper bound to match the specified
     reference image grid. This operation ignores differences in image
     transformation between input and reference image.

  -uniform number
     pad or crop the input image by a uniform number of voxels on all sides

  -mask image
     crop the input image according to the spatial extent of a mask image. The
     mask must share a common voxel grid with the input image but differences
     in image transformations are ignored. Note that even though only 3
     dimensions are cropped when using a mask, the bounds are computed by
     checking the extent for all dimensions. Note that by default a gap of 1
     voxel is left at all edges of the image to allow valid trilinear
     interpolation. This gap can be modified with the -uniform option but by
     default it does not extend beyond the FOV unless -crop_unbound is used.

  -crop_unbound
     Allow padding beyond the original FOV when cropping.

  -axis index spec  (multiple uses permitted)
     pad or crop the input image along the provided axis (defined by index).
     The specifier argument defines the number of voxels added or removed on
     the lower or upper end of the axis (-axis index delta_lower,delta_upper)
     or acts as a voxel selection range (-axis index start:stop). In both
     modes, values are relative to the input image (overriding all other
     extent-specifying options). Negative delta specifier values trigger the
     inverse operation (pad instead of crop and vice versa) and negative range
     specifier trigger padding. Note that the deprecated commands 'mrcrop' and
     'mrpad' used range-based and delta-based -axis indices, respectively.

  -all_axes
     Crop or pad all, not just spatial axes.

General options

  -fill number
     Use number as the out of bounds value. nan, inf and -inf are valid
     arguments. (Default: 0.0)

Stride options

  -strides spec
     specify the strides of the output data in memory; either as a
     comma-separated list of (signed) integers, or as a template image from
     which the strides shall be extracted and used. The actual strides produced
     will depend on whether the output image format can support it.

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.

解像度の変更する場合の基本的な使い方は、以下の通り。

mrgrid <入力画像> regrid -voxel <値> <出力画像>  # ボクセルサイズを指定
mrgrid <入力画像> regrid -scale <値> <出力画像>  # スケールを指定
mrgrid <入力画像> regrid -template <目的の解像度を持つ画像> <出力画像>  # 目的の解像度を持つ画像を指定

3. 使用例

3D-T1WI(T1w.nii.gz)の解像度を変更する。

3D-T1WI(T1w.nii.gz)の解像度を確認してみる。

mrinfo T1w.nii.gz
************************************************
Image name:          "T1w.nii.gz"
************************************************
  Dimensions:        192 x 256 x 256
  Voxel size:        0.9 x 0.9375 x 0.9375
  Data strides:      [ -1 2 3 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         signed 16 bit integer (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:               0.9998     0.01794   0.0003439      -82.89
                         -0.01788      0.9946      0.1023      -113.6
                         0.001492     -0.1023      0.9948      -114.6
  comments:          6.0.3:b862cdd5

3.1. ボクセルサイズを指定(オプション:-voxel

-voxelオプションを用いて、以下のコマンドを実行。

ボクセルサイズを1mm isotropicにする。

mrgrid T1w.nii.gz regrid -voxel 1 T1w_1mm_iso.nii.gz

解像度を確認してみる。ボクセルサイズが1 x 1 x 1(1mm iso)になっている

mrinfo T1w_1mm_iso.nii.gz
************************************************
Image name:          "T1w_1mm_iso.nii.gz"
************************************************
  Dimensions:        173 x 240 x 240
  Voxel size:        1 x 1 x 1
  Data strides:      [ -1 2 3 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:               0.9998     0.01794   0.0003439      -82.94
                         -0.01788      0.9946      0.1023      -113.6
                         0.001492     -0.1023      0.9948      -114.5
  comments:          6.0.3:b862cdd5
  mrtrix_version:    3.0.0-40-g3e1ed225

3.2. スケールを指定(オプション:-scale

-scaleオプションを用いて、以下のコマンドを実行。

スケールを2にして、解像度を2倍にする。

mrgrid T1w.nii.gz regrid -scale 2 T1w_scale2.nii.gz

解像度を確認してみる。解像度が173 x 240 x 240からになっている。

mrinfo T1w_scale2.nii.gz
************************************************
Image name:          "T1w_scale2.nii.gz"
************************************************
  Dimensions:        384 x 512 x 512
  Voxel size:        0.45 x 0.46875 x 0.46875
  Data strides:      [ -1 2 3 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:               0.9998     0.01794   0.0003439      -83.12
                         -0.01788      0.9946      0.1023      -113.9
                         0.001492     -0.1023      0.9948      -114.8
  comments:          6.0.3:b862cdd5
  mrtrix_version:    3.0.0-40-g3e1ed225

3.3. ボクセルサイズを指定(オプション:-voxel

-voxelオプションを用いて、以下のコマンドを実行。

ボクセルサイズを1mm isotropicにする。

mrgrid T1w.nii.gz regrid -voxel 1 T1w_1mm_iso.nii.gz

解像度を確認してみる。ボクセルサイズが1 x 1 x 1(1mm iso)になっている。

mrinfo T1w_1mm_iso.nii.gz
************************************************
Image name:          "T1w_1mm_iso.nii.gz"
************************************************
  Dimensions:        173 x 240 x 240
  Voxel size:        1 x 1 x 1
  Data strides:      [ -1 2 3 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:               0.9998     0.01794   0.0003439      -82.94
                         -0.01788      0.9946      0.1023      -113.6
                         0.001492     -0.1023      0.9948      -114.5
  comments:          6.0.3:b862cdd5
  mrtrix_version:    3.0.0-40-g3e1ed225

3.4. 目的の解像度を持つ画像を指定(オプション:-template

標準脳(MNI152)の3D-T1WI(MNI152_T1_2mm.nii.gz)と同じ解像度にする。標準脳(MNI152_T1_2mm.nii.gz)の解像度は以下。

mrinfo MNI152_T1_2mm.nii.gz
************************************************
Image name:          "MNI152_T1_2mm.nii.gz"
************************************************
  Dimensions:        91 x 109 x 91
  Voxel size:        2 x 2 x 2
  Data strides:      [ -1 2 3 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         signed 16 bit integer (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:                    1           0           0         -90
                               -0           1           0        -126
                               -0           0           1         -72
  comments:          FSL5.0

個人脳(T1w.nii.gz)を標準脳(MNI152_T1_2mm.nii.gz)の解像度に合わせるには、-templateオプションを用いて、以下のコマンドを実行。

mrgrid T1w.nii.gz regrid -template MNI152_T1_2mm.nii.gz T1w_MNIreso.nii.gz

解像度を確認してみる。解像度が標準脳(MNI152_T1_2mm.nii.gz)と同じになっている。

mrinfo T1w_MNIreso.nii.gz
************************************************
Image name:          "T1w_MNIreso.nii.gz"
************************************************
  Dimensions:        91 x 109 x 91
  Voxel size:        2 x 2 x 2
  Data strides:      [ -1 2 3 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:                    1           0           0         -90
                               -0           1           0        -126
                               -0           0           1         -72
  comments:          6.0.3:b862cdd5
  mrtrix_version:    3.0.0-40-g3e1ed225

【FSL/MRtrix】4D画像から3D画像を抽出


1. 目的
2. FSLを用いる場合
2.1. コマンド
2.2. 使用例
3. MRtrixを用いる場合
3.1. コマンド
3.2. 使用例


1. 目的

  • 4D画像から3D画像を抽出

2. FSLを用いる場合

2.1. コマンド

FSLfslroiコマンドを用いる。

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

クリックして展開
Usage: fslroi <input> <output> <xmin> <xsize> <ymin> <ysize> <zmin> <zsize>
       fslroi <input> <output> <tmin> <tsize>

       fslroi <input> <output> <xmin> <xsize> <ymin> <ysize> <zmin> <zsize> <tmin> <tsize>
Note: indexing (in both time and space) starts with 0 not 1! Inputting -1 for a size will set it to the full image extent for that dimension.

4D画像から3D画像を抽出する際の、基本的な使い方は以下の通り。

fslroi <入力画像> <出力画像> <Volume Index> <Volume Indexから残したいVolume数>

2.2. 使用例

例えば、5ttgen等で作成した以下のような5つの組織画像(5tt.nii.gz)が4D画像となっている場合。

Pathological tissue(Volume 4th)を取り除くには、次のようにコマンドを実行する。FSLではVolumeのIndexを0から数える。つまり、1番目のVolumeのIndexは0となる。以下のコードを翻訳すると、「Volume Index0番から数えて4 Volumesまでを残す」ということになる。

fslroi 5tt.nii.gz 4tt.nii.gz 0 4

fslhdコマンドを用いて、ボリューム数を確認すると、処理前で5 Volumesだったのが処理後に4 Volumesになっていることが分かる。使い方の詳細は、こちらの記事を参考に。

fslhd 5tt.nii.gz |grep ^dim4
fslhd 4tt.nii.gz |grep ^dim4
dim4        5  # 5tt.nii.gz
dim4        4  # 4tt.nii.gz

3. MRtrixを用いる場合

3.1. コマンド

MRtrixmrconvertコマンドを用いる。

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

クリックして展開
SYNOPSIS

     Perform conversion between different file types and optionally extract a
     subset of the input image

USAGE

     mrconvert [ options ] input output

        input        the input image.

        output       the output image.


DESCRIPTION

     If used correctly, this program can be a very useful workhorse. In
     addition to converting images between different formats, it can be used to
     extract specific studies from a data set, extract a specific region of
     interest, or flip the images. Some of the possible operations are
     described in more detail below.

     Note that for both the -coord and -axes options, indexing starts from 0
     rather than 1. E.g. -coord 3 <#> selects volumes (the fourth dimension)
     from the series; -axes 0,1,2 includes only the three spatial axes in the
     output image.

     Additionally, for the second input to the -coord option and the -axes
     option, you can use any valid number sequence in the selection, as well as
     the 'end' keyword (see the main documentation for details); this can be
     particularly useful to select multiple coordinates.

     The -vox option is used to change the size of the voxels in the output
     image as reported in the image header; note however that this does not
     re-sample the image based on a new voxel size (that is done using the
     mrresize command).

     By default, the intensity scaling parameters in the input image header are
     passed through to the output image header when writing to an integer
     image, and reset to 0,1 (i.e. no scaling) for floating-point and binary
     images. Note that the -scaling option will therefore have no effect for
     floating-point or binary output images.

     The -axes option specifies which axes from the input image will be used to
     form the output image. This allows the permutation, omission, or addition
     of axes into the output image. The axes should be supplied as a
     comma-separated list of axis indices. If an axis from the input image is
     to be omitted from the output image, it must either already have a size of
     1, or a single coordinate along that axis must be selected by the user by
     using the -coord option. Examples are provided further below.

     The -bvalue_scaling option controls an aspect of the import of diffusion
     gradient tables. When the input diffusion-weighting direction vectors have
     norms that differ substantially from unity, the b-values will be scaled by
     the square of their corresponding vector norm (this is how multi-shell
     acquisitions are frequently achieved on scanner platforms). However in
     some rare instances, the b-values may be correct, despite the vectors not
     being of unit norm (or conversely, the b-values may need to be rescaled
     even though the vectors are close to unit norm). This option allows the
     user to control this operation and override MRrtix3's automatic detection.

EXAMPLE USAGES

     Extract the first volume from a 4D image, and make the output a 3D image:
       $ mrconvert in.mif -coord 3 0 -axes 0,1,2 out.mif
     The -coord 3 0 option extracts, from axis number 3 (which is the fourth
     axis since counting begins from 0; this is the axis that steps across
     image volumes), only coordinate number 0 (i.e. the first volume). The
     -axes 0,1,2 ensures that only the first three axes (i.e. the spatial axes)
     are retained; if this option were not used in this example, then image
     out.mif would be a 4D image, but it would only consist of a single volume,
     and mrinfo would report its size along the fourth axis as 1.

     Extract slice number 24 along the AP direction:
       $ mrconvert volume.mif slice.mif -coord 1 24
     MRtrix3 uses a RAS (Right-Anterior-Superior) axis convention, and
     internally reorients images upon loading in order to conform to this as
     far as possible. So for non-exotic data, axis 1 should correspond
     (approximately) to the anterior-posterior direction.

     Extract only every other volume from a 4D image:
       $ mrconvert all.mif every_other.mif -coord 3 1:2:end
     This example demonstrates two features: Use of the colon syntax to
     conveniently specify a number sequence (in the format 'start:step:stop');
     and use of the 'end' keyword to generate this sequence up to the size of
     the input image along that axis (i.e. the number of volumes).

     Alter the image header to report a new isotropic voxel size:
       $ mrconvert in.mif isotropic.mif -vox 1.25
     By providing a single value to the -vox option only, the specified value
     is used to set the voxel size in mm for all three spatial axes in the
     output image.

     Alter the image header to report a new anisotropic voxel size:
       $ mrconvert in.mif anisotropic.mif -vox 1,,3.5
     This example will change the reported voxel size along the first and third
     axes (ideally left-right and inferior-superior) to 1.0mm and 3.5mm
     respectively, and leave the voxel size along the second axis (ideally
     anterior-posterior) unchanged.

     Turn a single-volume 4D image into a 3D image:
       $ mrconvert 4D.mif 3D.mif -axes 0,1,2
     Sometimes in the process of extracting or calculating a single 3D volume
     from a 4D image series, the size of the image reported by mrinfo will be
     "X x Y x Z x 1", indicating that the resulting image is in fact also 4D,
     it just happens to contain only one volume. This example demonstrates how
     to convert this into a genuine 3D image (i.e. mrinfo will report the size
     as "X x Y x Z".

     Insert an axis of size 1 into the image:
       $ mrconvert XYZD.mif XYZ1D.mif -axes 0,1,2,-1,3
     This example uses the value -1 as a flag to indicate to mrconvert where a
     new axis of unity size is to be inserted. In this particular example, the
     input image has four axes: the spatial axes X, Y and Z, and some form of
     data D is stored across the fourth axis (i.e. volumes). Due to insertion
     of a new axis, the output image is 5D: the three spatial axes (XYZ), a
     single volume (the size of the output image along the fourth axis will be
     1), and data D will be stored as volume groups along the fifth axis of the
     image.

     Manually reset the data scaling parameters stored within the image header
     to defaults:
       $ mrconvert with_scaling.mif without_scaling.mif -scaling 0.0,1.0
     This command-line option alters the parameters stored within the image
     header that provide a linear mapping from raw intensity values stored in
     the image data to some other scale. Where the raw data stored in a
     particular voxel is I, the value within that voxel is interpreted as:
     value = offset + (scale x I).  To adjust this scaling, the relevant
     parameters must be provided as a comma-separated 2-vector of
     floating-point values, in the format "offset,scale" (no quotation marks).
     This particular example sets the offset to zero and the scale to one,
     which equates to no rescaling of the raw intensity data.

Options for manipulating fundamental image properties

  -coord axis selection  (multiple uses permitted)
     retain data from the input image only at the coordinates specified in the
     selection along the specified axis. The selection argument expects a
     number sequence, which can also include the 'end' keyword.

  -vox sizes
     change the voxel dimensions reported in the output image header

  -axes axes
     specify the axes from the input image that will be used to form the output
     image

  -scaling values
     specify the data scaling parameters used to rescale the intensity values

Options for handling JSON (JavaScript Object Notation) files

  -json_import file
     import data from a JSON file into header key-value pairs

  -json_export file
     export data from an image header key-value pairs into a JSON file

Options to modify generic header entries

  -clear_property key  (multiple uses permitted)
     remove the specified key from the image header altogether.

  -set_property key value  (multiple uses permitted)
     set the value of the specified key in the image header.

  -append_property key value  (multiple uses permitted)
     append the given value to the specified key in the image header (this adds
     the value specified as a new line in the header value).

  -copy_properties source
     clear all generic properties and replace with the properties from the
     image / file specified.

Stride options

  -strides spec
     specify the strides of the output data in memory; either as a
     comma-separated list of (signed) integers, or as a template image from
     which the strides shall be extracted and used. The actual strides produced
     will depend on whether the output image format can support it.

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.

DW gradient table import options

  -grad file
     Provide the diffusion-weighted gradient scheme used in the acquisition in
     a text file. This should be supplied as a 4xN text file with each line is
     in the format [ X Y Z b ], where [ X Y Z ] describe the direction of the
     applied gradient, and b gives the b-value in units of s/mm^2. If a
     diffusion gradient scheme is present in the input image header, the data
     provided with this option will be instead used.

  -fslgrad bvecs bvals
     Provide the diffusion-weighted gradient scheme used in the acquisition in
     FSL bvecs/bvals format files. If a diffusion gradient scheme is present in
     the input image header, the data provided with this option will be instead
     used.

  -bvalue_scaling mode
     enable or disable scaling of diffusion b-values by the square of the
     corresponding DW gradient norm (see Desciption). Valid choices are yes/no,
     true/false, 0/1 (default: automatic).

DW gradient table export options

  -export_grad_mrtrix path
     export the diffusion-weighted gradient table to file in MRtrix format

  -export_grad_fsl bvecs_path bvals_path
     export the diffusion-weighted gradient table to files in FSL (bvecs /
     bvals) format

Options for importing phase-encode tables

  -import_pe_table file
     import a phase-encoding table from file

  -import_pe_eddy config indices
     import phase-encoding information from an EDDY-style config / index file
     pair

Options for exporting phase-encode tables

  -export_pe_table file
     export phase-encoding table to file

  -export_pe_eddy config indices
     export phase-encoding information to an EDDY-style config / index file
     pair

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.

4D画像から3D画像を抽出する際の、基本的な使い方は以下の通り。

mrconvert <入力画像> <出力画像> -coord <軸番号> <残したいボリューム数>

3.2. 使用例

例えば、5ttgen等で作成した以下のような5つの組織画像(5tt.nii.gz)が4D画像となっている場合。

Pathological tissue(Volume 4th)を取り除くには、次のようにコマンドを実行する。MRtrixでもFSLと同様に、VolumeのIndexを0から数える。つまり、1番目のVolumeのIndexは0となる。また軸番号は、x, y, z, tの順番に0, 1, 2, 3であり、Volume数を操作するには、t軸(-coord 3)を操作することになる。以下のコードを翻訳すると、「Volume Index0番からVolume Index3番までを残す」ということになる。

mrconvert 5tt.nii.gz 4tt.nii.gz -coord 3 0:3

mrinfoコマンドを用いて、ボリューム数を確認すると、処理前で5 Volumesだったのが処理後に4 Volumesになっていることが分かる。使い方の詳細は、こちらの記事を参考に。

mrinfo 5tt.nii.gz 4tt.nii.gz
************************************************
Image name:          "5tt.nii.gz"
************************************************
  Dimensions:        168 x 185 x 169 x 5
  Voxel size:        0.9 x 0.9375 x 0.9375 x ?
  Data strides:      [ 1 2 3 4 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:               0.9998     0.01794   0.0003439      -70.81
                         -0.01788      0.9946      0.1023       -88.1
                         0.001492     -0.1023      0.9948      -56.89
  comments:          6.0.3:b862cdd5
  mrtrix_version:    3.0.0-40-g3e1ed225
************************************************
Image name:          "4tt.nii.gz"
************************************************
  Dimensions:        168 x 185 x 169 x 4
  Voxel size:        0.9 x 0.9375 x 0.9375 x ?
  Data strides:      [ 1 2 3 4 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:               0.9998     0.01794   0.0003439      -70.81
                         -0.01788      0.9946      0.1023       -88.1
                         0.001492     -0.1023      0.9948      -56.89
  comments:          6.0.3:b862cdd5
  mrtrix_version:    3.0.0-40-g3e1ed225

【FSL/MRtrix】画像の切り取り・マスキング ~Masking~

目的

  • 画像の切り取り・マスキング ~Masking~

FSLを用いる場合

コマンド

FSLで画像の切り取り・マスキングをするには、fslmaths-masオプションを使用する。

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 <入力画像> -mas <マスク画像> <出力画像>

使用例

頭蓋除去されていないFA画像とマスク画像(緑)を、重ね合わせて表示した画像を以下に示す。

頭蓋除去されていないFA画像(FA.nii.gz)をマスク画像(mask.nii.gz)でマスキングするには、以下のコマンドを実行する。

fslmaths FA.nii.gz -mas mask.nii.gz FA_masked.nii.gz

マスキングして、頭蓋除去したFA画像は以下。

MRtrixを用いる場合

コマンド

MRtrixで画像の切り取り・マスキングをするには、mrcalc-multオプションを使用する。

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

クリックして展開
SYNOPSIS

     Apply generic voxel-wise mathematical operations to images

USAGE

     mrcalc [ options ] operand [ operand ... ]

        operand      an input image, intensity value, or the special keywords
                     'rand' (random number between 0 and 1) or 'randn' (random
                     number from unit std.dev. normal distribution) or the
                     mathematical constants 'e' and 'pi'.


DESCRIPTION

     This command will only compute per-voxel operations. Use 'mrmath' to
     compute summary statistics across images or along image axes.

     This command uses a stack-based syntax, with operators (specified using
     options) operating on the top-most entries (i.e. images or values) in the
     stack. Operands (values or images) are pushed onto the stack in the order
     they appear (as arguments) on the command-line, and operators (specified
     as options) operate on and consume the top-most entries in the stack, and
     push their output as a new entry on the stack.

     As an additional feature, this command will allow images with different
     dimensions to be processed, provided they satisfy the following
     conditions: for each axis, the dimensions match if they are the same size,
     or one of them has size one. In the latter case, the entire image will be
     replicated along that axis. This allows for example a 4D image of size [ X
     Y Z N ] to be added to a 3D image of size [ X Y Z ], as if it consisted of
     N copies of the 3D image along the 4th axis (the missing dimension is
     assumed to have size 1). Another example would a single-voxel 4D image of
     size [ 1 1 1 N ], multiplied by a 3D image of size [ X Y Z ], which would
     allow the creation of a 4D image where each volume consists of the 3D
     image scaled by the corresponding value for that volume in the
     single-voxel image.

EXAMPLE USAGES

     Double the value stored in every voxel:
       $ mrcalc a.mif 2 -mult r.mif
     This performs the operation: r = 2*a  for every voxel a,r in images a.mif
     and r.mif respectively.

     A more complex example:
       $ mrcalc a.mif -neg b.mif -div -exp 9.3 -mult r.mif
     This performs the operation: r = 9.3*exp(-a/b)

     Another complex example:
       $ mrcalc a.mif b.mif -add c.mif d.mif -mult 4.2 -add -div r.mif
     This performs: r = (a+b)/(c*d+4.2).

     Rescale the densities in a SH l=0 image:
       $ mrcalc ODF_CSF.mif 4 pi -mult -sqrt -div ODF_CSF_scaled.mif
     This applies the spherical harmonic basis scaling factor: 1.0/sqrt(4*pi),
     such that a single-tissue voxel containing the same intensities as the
     response function of that tissue should contain the value 1.0.

basic operations

  -abs  (multiple uses permitted)
     |%1| : return absolute value (magnitude) of real or complex number

  -neg  (multiple uses permitted)
     -%1 : negative value

  -add  (multiple uses permitted)
     (%1 + %2) : add values

  -subtract  (multiple uses permitted)
     (%1 - %2) : subtract nth operand from (n-1)th

  -multiply  (multiple uses permitted)
     (%1 * %2) : multiply values

  -divide  (multiple uses permitted)
     (%1 / %2) : divide (n-1)th operand by nth

  -min  (multiple uses permitted)
     min (%1, %2) : smallest of last two operands

  -max  (multiple uses permitted)
     max (%1, %2) : greatest of last two operands

comparison operators

  -lt  (multiple uses permitted)
     (%1 < %2) : less-than operator (true=1, false=0)

  -gt  (multiple uses permitted)
     (%1 > %2) : greater-than operator (true=1, false=0)

  -le  (multiple uses permitted)
     (%1 <= %2) : less-than-or-equal-to operator (true=1, false=0)

  -ge  (multiple uses permitted)
     (%1 >= %2) : greater-than-or-equal-to operator (true=1, false=0)

  -eq  (multiple uses permitted)
     (%1 == %2) : equal-to operator (true=1, false=0)

  -neq  (multiple uses permitted)
     (%1 != %2) : not-equal-to operator (true=1, false=0)

conditional operators

  -if  (multiple uses permitted)
     (%1 ? %2 : %3) : if first operand is true (non-zero), return second
     operand, otherwise return third operand

  -replace  (multiple uses permitted)
     (%1, %2 -> %3) : Wherever first operand is equal to the second operand,
     replace with third operand

power functions

  -sqrt  (multiple uses permitted)
     sqrt (%1) : square root

  -pow  (multiple uses permitted)
     %1^%2 : raise (n-1)th operand to nth power

nearest integer operations

  -round  (multiple uses permitted)
     round (%1) : round to nearest integer

  -ceil  (multiple uses permitted)
     ceil (%1) : round up to nearest integer

  -floor  (multiple uses permitted)
     floor (%1) : round down to nearest integer

logical operators

  -not  (multiple uses permitted)
     !%1 : NOT operator: true (1) if operand is false (i.e. zero)

  -and  (multiple uses permitted)
     (%1 && %2) : AND operator: true (1) if both operands are true (i.e.
     non-zero)

  -or  (multiple uses permitted)
     (%1 || %2) : OR operator: true (1) if either operand is true (i.e.
     non-zero)

  -xor  (multiple uses permitted)
     (%1 ^^ %2) : XOR operator: true (1) if only one of the operands is true
     (i.e. non-zero)

classification functions

  -isnan  (multiple uses permitted)
     isnan (%1) : true (1) if operand is not-a-number (NaN)

  -isinf  (multiple uses permitted)
     isinf (%1) : true (1) if operand is infinite (Inf)

  -finite  (multiple uses permitted)
     finite (%1) : true (1) if operand is finite (i.e. not NaN or Inf)

complex numbers

  -complex  (multiple uses permitted)
     (%1 + %2 i) : create complex number using the last two operands as
     real,imaginary components

  -polar  (multiple uses permitted)
     (%1 /_ %2) : create complex number using the last two operands as
     magnitude,phase components (phase in radians)

  -real  (multiple uses permitted)
     real (%1) : real part of complex number

  -imag  (multiple uses permitted)
     imag (%1) : imaginary part of complex number

  -phase  (multiple uses permitted)
     phase (%1) : phase of complex number (use -abs for magnitude)

  -conj  (multiple uses permitted)
     conj (%1) : complex conjugate

  -proj  (multiple uses permitted)
     proj (%1) : projection onto the Riemann sphere

exponential functions

  -exp  (multiple uses permitted)
     exp (%1) : exponential function

  -log  (multiple uses permitted)
     log (%1) : natural logarithm

  -log10  (multiple uses permitted)
     log10 (%1) : common logarithm

trigonometric functions

  -cos  (multiple uses permitted)
     cos (%1) : cosine

  -sin  (multiple uses permitted)
     sin (%1) : sine

  -tan  (multiple uses permitted)
     tan (%1) : tangent

  -acos  (multiple uses permitted)
     acos (%1) : inverse cosine

  -asin  (multiple uses permitted)
     asin (%1) : inverse sine

  -atan  (multiple uses permitted)
     atan (%1) : inverse tangent

hyperbolic functions

  -cosh  (multiple uses permitted)
     cosh (%1) : hyperbolic cosine

  -sinh  (multiple uses permitted)
     sinh (%1) : hyperbolic sine

  -tanh  (multiple uses permitted)
     tanh (%1) : hyperbolic tangent

  -acosh  (multiple uses permitted)
     acosh (%1) : inverse hyperbolic cosine

  -asinh  (multiple uses permitted)
     asinh (%1) : inverse hyperbolic sine

  -atanh  (multiple uses permitted)
     atanh (%1) : inverse hyperbolic tangent

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.

基本的な使い方は、以下の通り。入力画像とバイナリーマスク画像(二値画像)を掛け算することで、マスキングをする。

mrcalc <入力画像> <バイナリーマスク画像> -mult <出力画像>

使用例

頭蓋除去されていないFA画像とマスク画像(緑)を、重ね合わせて表示した画像を以下に示す。

頭蓋除去されていないFA画像(FA.nii.gz)をマスク画像(mask.nii.gz)でマスキングするには、以下のコマンドを実行する。

mrcalc FA.nii.gz  mask.nii.gz -mult FA_masked.nii.gz

マスキングして、頭蓋除去したFA画像は以下。

【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

【MRIcron/MRIcroGL】MRIcron/MRIcroGLを用いたノイズ除去とマスク画像の作成


1. 目的
2. MRIcronを用いる場合
3. MRIcroGLを用いる場合


1. 目的

  • MRIcron/MRIcroGLを用いたバイナリーマスク画像の作成

ここでは、拡散MRI(b=0, SE-EPI)の脳周囲にあるノイズ除去するため、マスク画像を生成しノイズを除去する方法を解説する。

まず、画像を見てみる。拡散MRI(b=0, SE-EPI)を信号値0-10のスケールで表示すると、以下のような画像が表示される。この内、脳周囲の白い点々(ごま塩ノイズ)は本来観測すべきではない信号、つまりノイズである。これから、このノイズを除去していく。

2. MRIcronを用いる場合

MRIcronの基本操作は、以下の記事を参考にするとよい。

脳のマスク画像を作るには、ツールタブの「Draw/Intensity filter」を選択。あるいは、「Ctrl + I」を押す。

脳実質が欠けないようにThresholdを設定する。設定ができたら「Save highlighted as NIfTI or VOI」を選択。

保存先を指定して保存する。この時NIfTI形式として保存しておくと、他の脳画像解析ソフトで扱いやすい。

保存した画像を開くと、Intensity filterでしきい値処理された画像が表示される。この時、画像はまだ二値化されていない状態である。

ここで、スケールを最初と同じ0-10に設定して、脳周囲のノイズを確認してみる。脳周囲のノイズを除去されていることが分かる。

この画像を二値化してバイナリーマスク画像を作成したい場合、ツールタブの「Draw/Advanced/Brain mask」を選択し、ファイル名を指定して保存する。

二値化されたバイナリーマスク画像が生成される。

3. MRIcroGLを用いる場合

MRIcroGLの基本操作は、以下の記事を参考にするとよい。

脳のマスク画像を作るには、ツールタブの「Draw/Advanced/Intensity Filter」を選択。

「Action: Add to Drawing」となっている状態で、脳実質が欠けないようにしきい値を設定し「Apply」をクリック。

Intensity Filterのしきい値処理が適用されて残った領域が、関心領域として設定される。

関心領域が設定されている状態で、上タブの「Draw/Advanced/Mask Image」から、「Delete/Preserve regions with VOI」を選択。

ここで、スケールを最初と同じ0-10に設定して、脳周囲のノイズを確認してみる。脳周囲のノイズを除去されていることが分かる。

この画像を二値化してバイナリーマスク画像を作成したい場合、ツールタブの「Draw/Save VOI」を選択し、ファイル名を指定して保存する(ショートカットキー:CTRL+S)。

二値化されたバイナリーマスク画像が生成される。

【FSL】 FSLを用いた定量値の計測 ~Sampling~


1. 目的
2. コマンド
3. 使用例
3.1. 最小値と最大値
3.2. ボクセル数および容積
3.3. 平均値と標準偏差
3.4. マスク画像を用いた計測


1. 目的

  • 定量値(容積や拡散定量値など)の算出

2. コマンド

FSLfslstatsコマンドを用いて、定量値を算出することが可能。

fslstatsのヘルプは次の通り。

Usage: fslstats [preoptions] <input> [options]

preoption -t will give a separate output line for each 3D volume of a 4D timeseries
preoption -K < indexMask > will generate seperate n submasks from indexMask, for indexvalues 1..n where n is the maximum index value in indexMask, and generate statistics for each submask
Note - options are applied in order, e.g. -M -l 10 -M will report the non-zero mean, apply a threshold and then report the new nonzero mean

-l <lthresh> : set lower threshold
-u <uthresh> : set upper threshold
-r           : output <robust min intensity> <robust max intensity>
-R           : output <min intensity> <max intensity>
-e           : output mean entropy ; mean(-i*ln(i))
-E           : output mean entropy (of nonzero voxels)
-v           : output <voxels> <volume>
-V           : output <voxels> <volume> (for nonzero voxels)
-m           : output mean
-M           : output mean (for nonzero voxels)
-s           : output standard deviation
-S           : output standard deviation (for nonzero voxels)
-w           : output smallest ROI <xmin> <xsize> <ymin> <ysize> <zmin> <zsize> <tmin> <tsize> containing nonzero voxels
-x           : output co-ordinates of maximum voxel
-X           : output co-ordinates of minimum voxel
-c           : output centre-of-gravity (cog) in mm coordinates
-C           : output centre-of-gravity (cog) in voxel coordinates
-p <n>       : output nth percentile (n between 0 and 100)
-P <n>       : output nth percentile (for nonzero voxels)
-a           : use absolute values of all image intensities
-n           : treat NaN or Inf as zero for subsequent stats
-k <mask>    : use the specified image (filename) for masking - overrides lower and upper thresholds
-d <image>   : take the difference between the base image and the image specified here
-h <nbins>   : output a histogram (for the thresholded/masked voxels only) with nbins
-H <nbins> <min> <max>   : output a histogram (for the thresholded/masked voxels only) with nbins and histogram limits of min and max

Note - thresholds are not inclusive ie lthresh<allowed<uthresh

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

# 基本
fslstats  <入力画像> [オプション]

# マスクを適用する場合
fslstats  <入力画像> -mas <マスク画像> [オプション]

3. 使用例

fslstatsコマンドで、よくある使用例を紹介する。

3.1. 最小値と最大値

最小値と最大値は、オプション-Rを用いて計測する。

以下では、頭蓋除去済みの脳画像(T1_skull_stripped.nii.gz)における最小値と最大値を計測している。

fslstats T1_skull_stripped.nii.gz -R
0.000000 1307.000000 

3.2. ボクセル数および容積

ボクセル数および容積の計測は、オプション-Vを用いる。似ているオプションで-vがあるが、オプション-Vでは、信号値が0あるいはnanでない領域を対象に、計測する。他にも小文字・大文字で区別しているオプションがあるが、信号値がある領域のみを対象にしたい場合は、基本的に大文字オプション(例:-V, -M, -S)で計測するとよい。

以下では、灰白質(GM_seg.nii.gz)の容積を計測している。

fslstats GM_seg.nii.gz -V
882175 697814.250000  # 左からボクセル数、容積(mm^3)

3.3. 平均値と標準偏差

平均値を計測するにはオプション-Mを、標準偏差を計測するにはオプション-Sを用いる。

以下では、FA(FA.nii.gz)の平均値を算出している。

fslstats FA.nii.gz -M
0.276669

以下では、FA(FA.nii.gz)の標準偏差を算出している。

fslstats FA.nii.gz -S
0.186857

3.4. マスク画像を用いた計測

計測したい領域がある場合、オプション-kで領域(マスク画像)を指定して計測するとよい。

例えば、白質(WM_seg.nii.gz)領域における、FAの平均値を計測したい場合、次のようになる。

fslstats FA.nii.gz -k WM_seg.nii.gz -M
0.507313

計測したい領域が複数あり、それらの領域にインデックス(値)が割り振られている画像(例:CSF=1, GM=2, WM=3の画像)があるとき、オプション-Kを用いると便利である。

例えば、白質を48領域分割したアトラス(JHU-ICBM-labels-1mm_indiv.nii.gz)を用いて平均値を計測する場合、次のようになる。

fslstats -K JHU-ICBM-labels-1mm_indiv.nii.gz FA.nii.gz  -M
0.491841 
0.512530 
0.591286 
0.667444 
... [省略]
0.294148

【FSL】FDT pipelineを用いた標準空間(MNI空間)への位置合わせ


1. 目的
2. 必要なファイル
3. 実行
4. 実際に実行されているコマンド
5. 出力画像


1. 目的

  • FDT pipelineを用いて、個人のDiffusion画像と標準空間(MNI空間)の構造画像の位置合わせ

2. 必要なファイル

必要なファイルは次の通り。

BEDPOSTXの使い方はこちら

.
├── DTI.bedpostX  # BEDPOSTXの出力フォルダ
├── T1.nii.gz  # BET前のBrain 3D-T1WI (FNIRT用)
└── T1_brain.nii.gz  # BET後のBrain 3D-T1WI (FLIRT用)

これに加えて、b0画像から脳を抽出した画像、nodif_brain.nii.gz を DTI.bedpostX にコピーしておく必要がある。

3. 実行

ターミナル(端末)でfslと入力。

fsl

Windowが立ち上がったら、「FDT diffusion」を選択。

各項目ごとにファイルを選択。注意すべきことは次の通り。

  • 「Main structural image」に指定する画像はBET後のBrain 3D-T1WI
  • 「Non-betted structural」に指定する画像は、BET前のBrain 3D-T1WI
  • Normal searchを「Full search」に変更

以上の設定ができたら、「Go」を選択して位置合わせを実行する。

4. 実際に実行されているコマンド

上記の操作を実行すると、ターミナル上にFSLのコマンドが生成され位置合わせが実行される。その時のコマンドは次の通り。

flirt -in DTI.bedpostX/nodif_brain \
    -ref T1_brain.nii.gz \
    -omat DTI.bedpostX/xfms/diff2str.mat \
    -searchrx -180 180 -searchry -180 180 -searchrz -180 180 \
    -dof 6 -cost corratio

convert_xfm -omat DTI.bedpostX/xfms/str2diff.mat \
    -inverse DTI.bedpostX/xfms/diff2str.mat

flirt -in T1_brain.nii.gz \
    -ref /opt/fsl/data/standard/MNI152_T1_2mm_brain \
    -omat DTI.bedpostX/xfms/str2standard.mat \
    -searchrx -180 180 -searchry -180 180 -searchrz -180 180 \
    -dof 12 -cost corratio

convert_xfm -omat DTI.bedpostX/xfms/standard2str.mat \
    -inverse DTI.bedpostX/xfms/str2standard.mat

convert_xfm -omat DTI.bedpostX/xfms/diff2standard.mat \
    -concat DTI.bedpostX/xfms/str2standard.mat DTI.bedpostX/xfms/diff2str.mat

convert_xfm -omat DTI.bedpostX/xfms/standard2diff.mat \
    -inverse DTI.bedpostX/xfms/diff2standard.mat

fnirt --in=T1.nii.gz \
    --aff=DTI.bedpostX/xfms/str2standard.mat \
    --cout=DTI.bedpostX/xfms/str2standard_warp \
    --config=T1_2_MNI152_2mm

invwarp -w DTI.bedpostX/xfms/str2standard_warp \
    -o DTI.bedpostX/xfms/standard2str_warp \
    -r T1_brain.nii.gz

convertwarp -o DTI.bedpostX/xfms/diff2standard_warp \
    -r /opt/fsl/data/standard/MNI152_T1_2mm \
    -m DTI.bedpostX/xfms/diff2str.mat \
    -w DTI.bedpostX/xfms/str2standard_warp

convertwarp -o DTI.bedpostX/xfms/standard2diff_warp \
    -r DTI.bedpostX/nodif_brain_mask \
    -w DTI.bedpostX/xfms/standard2str_warp \
    --postmat=DTI.bedpostX/xfms/str2diff.mat

5. 出力画像

処理が終わると、BEDPOSTXの出力フォルダ(DTI.bedpostX)に位置合わせの出力ファイルが保存される。

DTI.bedpostX/xfms/
├── diff2standard.mat
├── diff2standard_warp.nii.gz
├── diff2str.mat
├── eye.mat
├── standard2diff.mat
├── standard2diff_warp.nii.gz
├── standard2str.mat
├── standard2str_warp.nii.gz
├── str2diff.mat
├── str2standard.mat
└── str2standard_warp.nii.gz

【FSL】BEDPOSTXの使い方


* 1. 目的
* 2. BEDPOSTX
* 3.


1. 目的

  • BEDPOSTXの利用方法の取得

2. BEDPOSTX

BEDPOSTXの実行には、次のようなファイルが必要。

さらに、ファイル名は次のようにしておく必要がある。

Sub001/
├── bvals  # DWIのGradient Table
├── bvecs  # DWIのGradient Table
├── data.nii.gz  # DWI
└── nodif_brain_mask.nii.gz  # b=0のマスク

BEDPOSTXは、次のコマンドで実行できる。

bedpostx Sub001
# Usage: bedpostx <subject directory> [options]

3.

次のように、データを用意する。

$ tree HC003/
HC003/
├── HC003.bval   # DWIのGradient Table
├── HC003.bvec  # DWIのGradient Table
├── drHC003.nii.gz  # DWI
└── maskdrHC003.nii.gz  # b=0のマスク

この時、BEDPOSTXを実行するために必要なファイルが揃っているかをbedpostx_datacheckで確認することができる。

$ bedpostx_datacheck HC003/
HC003//data does not exist
HC003//nodif_brain_mask does not exist
 num lines in HC003//bvals 
cat: HC003//bvals: No such file or directory
0
 num words in HC003//bvals 
cat: HC003//bvals: No such file or directory
0
 num lines in HC003//bvecs 
cat: HC003//bvecs: No such file or directory
0
 num words in HC003//bvecs 
cat: HC003//bvecs: No such file or directory
0

ファイル名を修正。

tree HC003
HC003/
├── bvals
├── bvecs
├── data.nii.gz
└── nodif_brain_mask.nii.gz

再度、bedpostx_datacheckを実行する。

$ bedpostx_datacheck HC003/
HC003//data
data_type	INT16
dim1		120
dim2		120
dim3		84
dim4		137
datatype	4
pixdim1		1.700000
pixdim2		1.700000
pixdim3		1.700000
pixdim4		0.000000
cal_max		0.000000
cal_min		0.000000
file_type	NIFTI-1+

HC003//nodif_brain_mask
data_type	INT16
dim1		120
dim2		120
dim3		84
dim4		1
datatype	4
pixdim1		1.700000
pixdim2		1.700000
pixdim3		1.700000
pixdim4		0.000000
cal_max		0.000000
cal_min		0.000000
file_type	NIFTI-1+

 num lines in HC003//bvals 
1
 num words in HC003//bvals 
137
 num lines in HC003//bvecs 
3
 num words in HC003//bvecs 
411

BEDPOSTXデータのチェックができたら、BEDPOSTXを実行する。

$ bedpostx HC003/
subjectdir is /home/neuro/Documents/Yuya_S/temp/bedpostx/bedpostx_dir/HC003
Making bedpostx directory structure
Queuing preprocessing stages
...

BEDPOSTXを実行すると、「.bedpostX」フォルダが生成されこれがBEDPOSTXの生成ファイルである。

$ ls
HC003  HC003.bedpostX

BEDPOSTXで出力されるファイルは次の通り。

  • merged_th<i>samples – 4D volume – Samples from the distribution on theta
  • merged_ph<i>samples – 4D volume – Samples from the distribution on phi
  • theta and phi together represent the principal diffusion direction in spherical polar co-ordinates
  • merged_f<i>samples – 4D volume – Samples from the distribution on anisotropic volume fraction (see technical report).
  • mean_th<i>samples – 3D Volume – Mean of distribution on theta
  • mean_ph<i>samples – 3D Volume – Mean of distribution on phi
  • mean_f<i>samples – 3D Volume – Mean of distribution on f anisotropy. Note that in each voxel, fibres are ordered according to a decreasing mean f-value
  • an_dsamples – 3D Volume – Mean of distribution on diffusivity d
  • mean_d_stdsamples – 3D Volume – Mean of distribution on diffusivity variance parameter d_std (not produced if –model=1)
  • mean_S0samples – 3D Volume – Mean of distribution on T2w baseline signal intensity S0
  • dyads<i> – Mean of PDD distribution in vector form. Note that this file can be loaded into FSLeyes for easy viewing of diffusion directions
  • dyads<i>_dispersion – 3D Volume – Uncertainty on the estimated fibre orientation. Characterizes how wide the orientation distribution is around the respective PDD.(how is this calculated?)
  • nodif_brain_mask – binary mask created from nodif_brain – copied from input directory

結果を確認するには、以下のコマンドを実行。

cd HC003.bedpostX
fsleyes mean_fsumsamples.nii.gz \
  dyads1.nii.gz         -ot linevector -xc 1 0 0 -yc 1 0 0 -zc 1 0 0 -lw 2 \
  dyads2_thr0.05.nii.gz -ot linevector -xc 0 1 0 -yc 0 1 0 -zc 0 1 0 -lw 2 \
  dyads3_thr0.05.nii.gz -ot linevector -xc 0 0 1 -yc 0 0 1 -zc 0 0 1 -lw 2

【PyTorch】サンプル⑨ 〜 動的グラフ 〜


1. 目的
2. 前準備
3. 計算グラフ
4. パッケージのインポート
5. モデルの定義
6. 使用するデータ
7. ニューラルネットワークのモデルの構築
8. 損失関数の定義
9. 最適化関数の定義
10. 学習回数
11. モデルへのデータ入出力 (Forward pass)
12. 損失(loss)の計算
13. 勾配の初期化
14. 勾配の計算
15. パラメータ(Weight)の更新
16. 実行
16.1. 9_dynamic_graph.py


1. 目的

  • PyTorchの特徴の一つである動的グラフに挑戦する。

2. 前準備

PyTorchのインストールはこちらから。

初めて、Google Colaboratoryを使いたい方は、こちらをご覧ください。

3. 計算グラフ

計算グラフには、静的グラフ動的グラフがあります。

TensorFlowやTheanoといったフレームワークは静的グラフを用いていますが、それに対してPyTorchは動的グラフを採用しています。

静的グラフと動的グラフの説明についてはこちらをご覧ください。

簡単に言うと、静的グラフは、事前にニューラルネットワーク構造を決めてから学習・推定を実行します。この時、一度決めたネットワーク構造は学習や推定時に変更することができません。

これに対し、動的グラフは、ネットワーク構造を定義する時点でネットワーク構造を固定する必要はなく、ネットワークに入力されるデータの様子をみて、ネットワークを構造を適宜変えることができます。故に、動的な計算グラフになります。

この動的グラフにより、柔軟な学習や推測が可能となります。

4. パッケージのインポート

import random
import torch

5. モデルの定義

【PyTorch】サンプル⑧ 〜 複雑なモデルの構築方法 〜で紹介した方法でニューラルネットワークモデルを構築していきます。

これまでは、レイヤ(層)を積み重ねるようにして、ニューラルネットワークモデルを定義していましたが、今回はネットワークモデルをクラスを使って定義します。クラスについては、こちらが参考になりました。

簡単のために、線形結合のみのネットワークを構築していきます。

クラスの名前をDynamicNetとします。引数torch.nn.Moduleはお決まりです。

class DynamicNet(torch.nn.Module):

コンストラクタで2層の線形結合層をインスタンス化します。

コンストラクタとインスタンスについては、こちらが参考になりました。

    def __init__(self, D_in, H, D_out):
        super(DynamicNet, self).__init__()
        self.input_linear = torch.nn.Linear(D_in, H)
        self.middle_linear = torch.nn.Linear(H, H)
        self.output_linear = torch.nn.Linear(H, D_out)

コンストラクタで定義された線形結合層をつかってモデルの予測値y_predを計算します。

テンソルの演算式は、forwardオブジェクトで自由に定義できます。

このチュートリアルでは、動的グラフの特徴を感じてもらうために、不思議なモデルを作っていきます。

    def forward(self, x):
        h_relu = self.input_linear(x).clamp(min=0)
        for _ in range(random.randint(0, 3)):
            h_relu = self.middle_linear(h_relu).clamp(min=0)
        y_pred = self.output_linear(h_relu)
        return y_pred

ポイントは、この部分です。

        for _ in range(random.randint(0, 3)):
            h_relu = self.middle_linear(h_relu).clamp(min=0)

この部分でやっていることは、forループの回数を0〜3回をランダムに決めて、中間層の線形結合層の数を随時変更しています。

6. 使用するデータ

バッチサイズNを64、入力の次元D_inを1000、隠れ層の次元Hを100、出力の次元D_outを10とします。

# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

入力xと予測したいyを乱数で定義します。

# Create random input and output data
x = torch.randn(N, D_in)
y = torch.randn(N, D_out)

7. ニューラルネットワークのモデルの構築

先程作成したニューラルネットワークモデルのクラスDynamicNetをつかってモデルmodelを定義します。

model = DynamicNet(D_in, H, D_out)

8. 損失関数の定義

二乗誤差をこのモデルの損失関数とします。

reductionのデフォルトはmeanですので、何も指定しなければtorch.nn.MSELossは平均二乗誤差を返します。

reduction=sumとした場合は、累積二乗誤差を算出します。

criterion = torch.nn.MSELoss(reduction='sum')

9. 最適化関数の定義

optimパッケージを使って最適化関数としてMomentum SGDを定義します。

lrは学習率です。

optimizer = torch.optim.SGD(model.parameters(), lr=1e-4, momentum=0.9)

10. 学習回数

学習回数は500回にします。

for t in range(500):

11. モデルへのデータ入出力 (Forward pass)

定義したニューラルネットワークモデルmodelへデータxを入力し、予測値y_predを取得します。

    y_pred = model(x)

12. 損失(loss)の計算

定義した損失関数Momentum SGDで予測値y_predと真値yとの間の損失を計算します。

    loss = criterion(y_pred, y)
    if t % 100 == 99:
        print(t, loss.item())

13. 勾配の初期化

逆伝播(backward)させる前に、モデルのパラメータが持つ勾配を0(ゼロ)で初期化します。

    optimizer.zero_grad()

14. 勾配の計算

backwardメソッドでモデルパラメータ(Weight)の勾配を算出します。

    loss.backward()

15. パラメータ(Weight)の更新

stepメソッドでモデルのパラメータを更新します。

    optimizer.step()

16. 実行

以下のコードを9_dynamic_graph.pyとして保存します。

16.1. 9_dynamic_graph.py

import random
import torch


class DynamicNet(torch.nn.Module):
    def __init__(self, D_in, H, D_out):
        super(DynamicNet, self).__init__()
        self.input_linear = torch.nn.Linear(D_in, H)
        self.middle_linear = torch.nn.Linear(H, H)
        self.output_linear = torch.nn.Linear(H, D_out)

    def forward(self, x):
        h_relu = self.input_linear(x).clamp(min=0)
        for _ in range(random.randint(0, 3)):
            h_relu = self.middle_linear(h_relu).clamp(min=0)
        y_pred = self.output_linear(h_relu)
        return y_pred


# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

# Create random Tensors to hold inputs and outputs
x = torch.randn(N, D_in)
y = torch.randn(N, D_out)

# Construct our model by instantiating the class defined above
model = DynamicNet(D_in, H, D_out)

# Construct our loss function and an Optimizer.
criterion = torch.nn.MSELoss(reduction='sum')
optimizer = torch.optim.SGD(model.parameters(), lr=1e-4, momentum=0.9)
for t in range(500):
    # Forward pass: Compute predicted y by passing x to the model
    y_pred = model(x)

    # Compute and print loss
    loss = criterion(y_pred, y)
    if t % 100 == 99:
        print(t, loss.item())

    # Zero gradients, perform a backward pass, and update the weights.
    optimizer.zero_grad()
    loss.backward()
    optimizer.step()

コードを保存したら、実行してみましょう。

左の数字が学習回数、右の数値がモデルの推定値と実際の答えと二乗誤差です。

python3 9_dynamic_graph.py
99 3.707427978515625
199 1.2958085536956787
299 1.837145209312439
399 0.9760080575942993
499 0.43565574288368225

参考のために、静的グラフの場合の結果を載せておきます。

動的グラフの方が損失の収束が早いことがわかります。

In:

python3 9_static_graph.py

Out:

99 2.5070412158966064
199 0.04455046355724335
299 0.0011508199386298656
399 3.597284376155585e-05
499 1.361027443635976e-06

【PyTorch】サンプル⑧ 〜 複雑なモデルの構築方法 〜


1. 目的
2. 前準備
3. PyTorchのインポート
4. モデルの定義
5. 使用するデータ
6. ニューラルネットワークのモデルを構築
7. 損失関数の定義
8. 最適化関数の定義
9. 学習回数
10. モデルへのデータ入出力 (Forward pass)
11. 損失(loss)の計算
12. 勾配の初期化
13. 勾配の計算
14. 実行
14.1. 8_class_model.py
15. 終わりに


1. 目的

このチュートリアルに至るまでは、ニューラルネットワークモデルの定義を積み木を積み重ねるように単純なシーケンスtorch.nn.Sequentialで構築していました。

このtorch.nn.Sequentialを用いた方法は、モデルの定義が簡単である反面、ネットワーク構造も簡素なものしか作ることができません。

例えば、torch.nn.Sequentialでは、ResNetネットワーク構造を構築することがきません。

(Credit:Deep Residual Learning for Image Recognition. Kaiming He, Xiangyu Zhang, Shaoqing Ren, Jian Sun. arXiv:1512.03385 [cs.CV]
(or arXiv:1512.03385v1 [cs.CV] for this version) )

このチュートリアルでは、PyTorch: Custom nn Modulesを参考に、より複雑なニューラルネットワークのモデルが構築できる方法を紹介します。

2. 前準備

PyTorchのインストールはこちらから。

初めて、Google Colaboratoryを使いたい方は、こちらをご覧ください。

3. PyTorchのインポート

import torch

4. モデルの定義

これまでは、レイヤ(層)を積み重ねるようにして、ニューラルネットワークモデルを定義していましたが、今回はネットワークモデルをクラスを使って定義します。クラスについては、こちらが参考になりました。

簡単のために、線形結合が2層のネットワークを構築していきます。

クラスの名前をTwoLayerNetとします。引数torch.nn.Moduleはお決まりです。

class TwoLayerNet(torch.nn.Module):

コンストラクタで2層の線形結合層をインスタンス化します。

コンストラクタとインスタンスについては、こちらが参考になりました。

    def __init__(self, D_in, H, D_out):
        super(TwoLayerNet, self).__init__()
        self.linear1 = torch.nn.Linear(D_in, H)
        self.linear2 = torch.nn.Linear(H, D_out)

コンストラクタで定義された線形結合層をつかってモデルの予測値y_predを計算します。

テンソルの演算式は、forwardオブジェクトで自由に定義できます。

    def forward(self, x):
        h_relu = self.linear1(x).clamp(min=0)
        y_pred = self.linear2(h_relu)
        return y_pred

5. 使用するデータ

バッチサイズNを64、入力の次元D_inを1000、隠れ層の次元Hを100、出力の次元D_outを10とします。

# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

入力xと予測したいyを乱数で定義します。

# Create random input and output data
x = torch.randn(N, D_in)
y = torch.randn(N, D_out)

6. ニューラルネットワークのモデルを構築

先程作成したニューラルネットワークモデルのクラスTwoLayerNetをつかってモデルを定義します。

model = TwoLayerNet(D_in, H, D_out)

7. 損失関数の定義

二乗誤差をこのモデルの損失関数とします。

reductionのデフォルトはmeanですので、何も指定しなければtorch.nn.MSELossは平均二乗誤差を返します。

reduction=sumとした場合は、累積二乗誤差を算出します。

criterion = torch.nn.MSELoss(reduction='sum')

8. 最適化関数の定義

optimパッケージを使って最適化関数として確率勾配降下法(SGD: stochastic gradient descent)を定義します。

lrは学習率です。

optimizer = torch.optim.SGD(model.parameters(), lr=1e-4)

9. 学習回数

学習回数は500回にします。

for t in range(500):

10. モデルへのデータ入出力 (Forward pass)

定義したニューラルネットワークモデルmodelへデータxを入力し、予測値y_predを取得します。

    y_pred = model(x)

11. 損失(loss)の計算

定義した損失関数SGDで予測値y_predと真値yとの間の損失を計算します。

    loss = criterion(y_pred, y)
    if t % 100 == 99:
        print(t, loss.item())

12. 勾配の初期化

逆伝播(backward)させる前に、モデルのパラメータが持つ勾配を0(ゼロ)で初期化します。

    optimizer.zero_grad()

13. 勾配の計算

backwardメソッドでモデルパラメータ(Weight)の勾配を算出します。

    loss.backward()

パラメータ(Weight)の更新

stepメソッドでモデルのパラメータを更新します。

    optimizer.step()

14. 実行

以下のコードを8_class_model.pyとして保存します。

14.1. 8_class_model.py

import torch


class TwoLayerNet(torch.nn.Module):
    def __init__(self, D_in, H, D_out):
        super(TwoLayerNet, self).__init__()
        self.linear1 = torch.nn.Linear(D_in, H)
        self.linear2 = torch.nn.Linear(H, D_out)

    def forward(self, x):
        h_relu = self.linear1(x).clamp(min=0)
        y_pred = self.linear2(h_relu)
        return y_pred


# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

# Create random Tensors to hold inputs and outputs
x = torch.randn(N, D_in)
y = torch.randn(N, D_out)

# Construct our model by instantiating the class defined above
model = TwoLayerNet(D_in, H, D_out)

# Construct our loss function and an Optimizer.
criterion = torch.nn.MSELoss(reduction='sum')
optimizer = torch.optim.SGD(model.parameters(), lr=1e-4)
for t in range(500):
    # Forward pass: Compute predicted y by passing x to the model
    y_pred = model(x)

    # Compute and print loss
    loss = criterion(y_pred, y)
    if t % 100 == 99:
        print(t, loss.item())

    # Zero gradients, perform a backward pass, and update the weights.
    optimizer.zero_grad()
    loss.backward()
    optimizer.step()

コードを保存したら、実行してみましょう。
In:

python3 8_class_model.py 

Out:

99 2.5070412158966064
199 0.04455046355724335
299 0.0011508199386298656
399 3.597284376155585e-05
499 1.361027443635976e-06

左の数字が学習回数、右の数値がパーセプトロンの推定値と実際の答えと二乗誤差です。

学習を重ねるごとに、二乗誤差が小さくなることがわかります。

15. 終わりに

今回は、クラスを使ってニューラルネットワークモデルを定義しました。

より複雑なモデルを構築する場合には、コンストラクタに線形結合層、畳み込み層、プーリング層などのコンポーネントを増やし、それらをforwardオブジェクトで、自由に組み替えてみてください。

【PyTorch】サンプル⑦ 〜 optim パッケージ 〜


1. 目的
2. 前準備
3. optim パッケージ
4. PyTorchのインポート
5. 使用するデータ
6. ニューラルネットワークのモデルを定義
7. 損失(loss)の定義
8. 学習パラメータ
9. 最適化関数(optimizer)の定義
10. 学習回数
11. モデルへのデータ入出力 (Forward pass)
12. 損失(loss)の計算
13. 勾配の初期化
14. 勾配の計算
15. パラメータ(Weight)の更新
16. 実行
16.1. 7_optim_package.py


1. 目的

  • PyTorch: optimを参考にPyTorchのoptimパッケージを使って最適化関数(optimizer)を定義する。

2. 前準備

PyTorchのインストールはこちらから。

初めて、Google Colaboratoryを使いたい方は、こちらをご覧ください。

3. optim パッケージ

optimは、最適化アルゴリズムの定義に用います。

これまでは、モデルのパラメータ(weight)を更新するために、自分の手で確率勾配降下法(SGD: stochastic gradient descent)のコードを作っていました。

今回のチュートリアルでは、optimパッケージを使ってパラメータを更新する最適化関数(optimizer)を定義していきます。

optimパッケージには、SGD+momentum、RMSProp, Adamなどディープラーニング界でよく使われる最適化アルゴリズムが複数あります。

4. PyTorchのインポート

import torch

5. 使用するデータ

バッチサイズNを64、入力の次元D_inを1000、隠れ層の次元Hを100、出力の次元D_outを10とします。

# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

入力(x)と予測したい(y)を乱数で定義します。

# Create random input and output data
x = torch.randn(N, D_in)
y = torch.randn(N, D_out)

6. ニューラルネットワークのモデルを定義

前回のnnパッケージのチュートリアルと同じように、ニューラルネットワークのモデルをnnパッケージを用いて定義します。

input > Linear(線型結合) > ReLU(活性化関数) > Linear(線型結合) > outputの順に層を積み重ねます。

model = torch.nn.Sequential(
    torch.nn.Linear(D_in, H),
    torch.nn.ReLU(),
    torch.nn.Linear(H, D_out),
)

7. 損失(loss)の定義

二乗誤差をnnパッケージを用いて計算します。

reductionのデフォルトはmeanですので、何も指定しなければtorch.nn.MSELossは平均二乗誤差を返します。

reduction=sumとした場合は、累積二乗誤差を算出します。

loss_fn = torch.nn.MSELoss(reduction='sum')

8. 学習パラメータ

学習率learning_rate1e-4とします。

learning_rate = 1e-4

9. 最適化関数(optimizer)の定義

最適化関数の定義は、torch.optimを使って簡単にできます。

このチュートリアルでは、最適化関数としてAdamを選択しています。

torch.optim.Adamの引数は、モデルのパラメータmodel.parameters()と学習率lr=learning_rateです。

optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate)

10. 学習回数

学習回数を500回とします。

for t in range(500):

11. モデルへのデータ入出力 (Forward pass)

定義したニューラルネットワークモデルmodelへデータxを入力し、予測値y_predを取得します。

    y_pred = model(x)

12. 損失(loss)の計算

定義した損失関数で予測値y_predと真値yとの間の損失を計算します。

    loss = loss_fn(y_pred, y)
    if t % 100 == 99:
        print(t, loss.item())

13. 勾配の初期化

逆伝播(backward)させる前に、モデルのパラメータが持つ勾配を0(ゼロ)で初期化します。

    optimizer.zero_grad()

(参考) optimパッケージを使わない場合、以下のように記述していました。

    model.zero_grad()

14. 勾配の計算

backwardメソッドでモデルパラメータ(Weight)の勾配を算出します。

    loss.backward()

15. パラメータ(Weight)の更新

stepメソッドでモデルのパラメータを更新します。

    optimizer.step()

(参考) optimパッケージを使わない場合、以下のように記述していました。

    with torch.no_grad():
        for param in model.parameters():
            param -= learning_rate param.grad

16. 実行

以下のコードを7_optim_package.pyとして保存します。

16.1. 7_optim_package.py

import torch

# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

# Create random Tensors to hold inputs and outputs
x = torch.randn(N, D_in)
y = torch.randn(N, D_out)

# Use the nn package to define our model and loss function.
model = torch.nn.Sequential(
    torch.nn.Linear(D_in, H),
    torch.nn.ReLU(),
    torch.nn.Linear(H, D_out),
)
loss_fn = torch.nn.MSELoss(reduction='sum')

# Use the optim package to define an Optimizer that will update the weights of
# the model for us.
learning_rate = 1e-4
optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate)
for t in range(500):
    # Forward pass: compute predicted y by passing x to the model.
    y_pred = model(x)

    # Compute and print loss.
    loss = loss_fn(y_pred, y)
    if t % 100 == 99:
        print(t, loss.item())

    # Before the backward pass, use the optimizer object to zero all of the
    # gradients for the variables it will update (which are the learnable
    # weights of the model).
    optimizer.zero_grad()

    # Backward pass: compute gradient of the loss with respect to model
    # parameters
    loss.backward()

    # Calling the step function on an Optimizer makes an update to its
    # parameters
    optimizer.step()

保存ができたら実行しましょう。

左の数字が学習回数、右の数値がパーセプトロンの推定値と実際の答えと二乗誤差です。

学習を重ねるごとに、二乗誤差が小さくなることがわかります。

In:

python3 7_optim_package.py 

Out:

99 54.53140640258789
199 0.6615686416625977
299 0.004912459757179022
399 3.133853169856593e-05
499 1.0647939063801459e-07

【PyTorch】サンプル⑥ 〜 nn パッケージ 〜


1. 目的
2. 前準備
3. nn パッケージ
4. PyTorchのインポート
5. 使用するデータ
6. ニューラルネットワークのモデルを定義
7. 損失(loss)の定義
8. 学習パラメータ
9. モデルへのデータ入出力 (Forward pass)
10. 損失(loss)の計算
11. 勾配の初期化
12. 勾配の計算
13. パラメータ(Weight)の更新
14. 実行
14.1. 6_nn_package.py
15. 終わりに


1. 目的

  • PyTorch: nnを参考にPyTorchのnnパッケージを扱う。
  • nnパッケージの便利さを感じる。

2. 前準備

PyTorchのインストールはこちらから。

初めて、Google Colaboratoryを使いたい方は、こちらをご覧ください。

3. nn パッケージ

nnは、ニューラルネットワークの構築に用いる。

PyTorchの自動微分autogradによって、計算グラフやパラメータの勾配を簡単に計算することができます。ですが、自動微分だけで複雑なニューラルネットワークを定義するのは困難です。そこで活躍するのがnnパッケージです。

nnパッケージを用いて、ニューラルネットワークを一つのモジュールとして定義することができます。

4. PyTorchのインポート

import torch

5. 使用するデータ

バッチサイズNを64、入力の次元D_inを1000、隠れ層の次元Hを100、出力の次元D_outを10とします。

# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

入力(x)と予測したい(y)を乱数で定義します。

# Create random input and output data
x = torch.randn(N, D_in)
y = torch.randn(N, D_out)

6. ニューラルネットワークのモデルを定義

ニューラルネットワークのモデルをnnパッケージを用いて定義します。

定義の仕方は、大きく2つありますがここでは一番簡単なtorch.nn.Sequentialを使います。

作り方は簡単で、任意の層を積み重ねていくだけです。この例では、input > Linear(線型結合) > ReLU(活性化関数) > Linear(線型結合) > outputの順に層が積み重なっています。

model = torch.nn.Sequential(
    torch.nn.Linear(D_in, H),
    torch.nn.ReLU(),
    torch.nn.Linear(H, D_out),
)

7. 損失(loss)の定義

二乗誤差もnnパッケージを用いて計算することができます。

reductionのデフォルトはmeanですので、何も指定しなければtorch.nn.MSELossは平均二乗誤差を返します。

reduction=sumとした場合は、累積二乗誤差を算出します。

loss_fn = torch.nn.MSELoss(reduction='sum')

(参考) nnパッケージを使う前は以下のように記述していた。

    loss = (y_pred - y).pow(2).sum()

8. 学習パラメータ

学習率を1e-4として、学習回数を500回とします。

learning_rate = 1e-4
for t in range(500):

9. モデルへのデータ入出力 (Forward pass)

定義したニューラルネットワークモデルへデータxを入力し、予測値y_predを取得します。

y_pred = model(x)

(参考) nnパッケージを使う前は以下のように記述していた。

    y_pred = x.mm(w1).clamp(min=0).mm(w2)

10. 損失(loss)の計算

定義した損失関数で予測値y_predと真値yとの間の損失を計算します。

    loss = loss_fn(y_pred, y)
    if t % 100 == 99:
        print(t, loss.item())

11. 勾配の初期化

逆伝播(backward)させる前に、モデルのパラメータが持つ勾配を0(ゼロ)で初期化します。

    model.zero_grad()

12. 勾配の計算

backwardメソッドでモデルパラメータ(Weight)の勾配を算出します。

    loss.backward()

13. パラメータ(Weight)の更新

確率勾配降下法(SGD: stochastic gradient descent)で、Weightを更新する。

    with torch.no_grad():
        for param in model.parameters():
            param -= learning_rate param.grad

(参考) nnパッケージを使う前は以下のように記述していた。

    with torch.no_grad():
        w1 -= learning_rate w1.grad
        w2 -= learning_rate w2.grad

14. 実行

以下のコードを6_nn_package.pyとして保存します。

14.1. 6_nn_package.py

import torch

# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

# Create random Tensors to hold inputs and outputs
x = torch.randn(N, D_in)
y = torch.randn(N, D_out)

# Use the nn package to define our model as a sequence of layers. 
model = torch.nn.Sequential(
    torch.nn.Linear(D_in, H),
    torch.nn.ReLU(),
    torch.nn.Linear(H, D_out),
)

# The nn package also contains definitions of popular loss functions; in this
# case we will use Mean Squared Error (MSE) as our loss function.
loss_fn = torch.nn.MSELoss(reduction='sum')

learning_rate = 1e-4
for t in range(500):
    # Forward pass: compute predicted y by passing x to the model. 
    y_pred = model(x)

    # Compute and print loss.
    loss = loss_fn(y_pred, y)
    if t % 100 == 99:
        print(t, loss.item())

    # Zero the gradients before running the backward pass.
    model.zero_grad()

    # Backward pass: compute gradient of the loss with respect to all the learnable
    loss.backward()

    # Update the weights using gradient descent.
    with torch.no_grad():
        for param in model.parameters():
            param -= learning_rate param.grad

保存ができたら実行しましょう。

左の数字が学習回数、右の数値がパーセプトロンの推定値と実際の答えと二乗誤差です。

学習を重ねるごとに、二乗誤差が小さくなることがわかります。

$ python3 6_nn_package.py 
99 2.5003600120544434
199 0.06977272033691406
299 0.003307548351585865
399 0.00018405442824587226
499 1.1299152902211063e-05

15. 終わりに

多層のニューラルネットワークには、膨大な量のパラメータが存在しています。

nnパッケージを用いる前は、各パラメータごとに勾配計算やパラメータの更新などをしていましたが、それでは記述が困難です。

nnパッケージを用いることで、楽に勾配計算やパラメータの更新等が実行できることが感じてもらえたら嬉しいです(^^)!。

【PyTorch】サンプル⑤ 〜Static Graphs(静的グラフ)〜


1. 目的
2. 前準備
3. 静的グラフ動的グラフ
4. パッケージのインポート
5. 使用するデータとプレースフォルダplaceholdersの用意
6. 重み付けの初期化
7. データ入力 (Forward pass)
8. 損失の計算
9. 重み(Weight)の勾配計算
10. 重み(w)の更新
11. 計算グラフの実行
12. 実行
12.1. 5_static_graphs.py


1. 目的

Deep Learning(ディープラーニング)で、よく使われるTensorFlowを使って計算グラフの設計からニューラルネットワークの学習をする。

PyTorchが採用している動的グラフを理解するために、静的グラフの代表であるTensorFlowに触れる。

2. 前準備

PyTorchのインストールはこちらから。

初めて、Google Colaboratoryを使いたい方は、こちらをご覧ください。

3. 静的グラフ動的グラフ

TensorFlowとPyTorchの大きな違いは、TensorFlowが静的計算グラフ(define-by-run)を採用しているのに対し、PyTorchは動的計算グラフ(define-and-run)を使用している点です。

静的グラフと動的グラフの説明についてはこちらをご覧ください。

簡潔に言うと、静的グラフは、事前にネットワーク構造を決めてから学習・推定を実行します。

これに対し、動的グラフは、ネットワーク構造を定義する時点でネットワーク構造を固定する必要はなく、ネットワークに入力されるデータの様子をみて、ネットワークを構造を適宜変えることができます。故に、動的な計算グラフになります。

この動的グラフにより、柔軟な学習や推測が可能となります。

4. パッケージのインポート

TensorFlowとNumPyをimportする。

import tensorflow as tf
import numpy as np

5. 使用するデータとプレースフォルダplaceholdersの用意

バッチサイズNを64、入力の次元D_inを1000、隠れ層の次元Hを100、出力の次元D_outを10とします。

# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

入力値xと予測したい値yの入れ物tf.placeholdersを用意する。

x = tf.placeholder(tf.float32, shape=(None, D_in))
y = tf.placeholder(tf.float32, shape=(None, D_out))

6. 重み付けの初期化

乱数random_normalメソッドを使って重みwを初期化します。

w1 = tf.Variable(tf.random_normal((D_in, H)))
w2 = tf.Variable(tf.random_normal((H, D_out)))

7. データ入力 (Forward pass)

入力xと重みw1を掛け算matmulすることで重み付けをしますh

重み付けした値(h)の要素から、maximum(h,tf.zeros(1))で、0以上のものは残し、0以下のものは0とします。

最後に、重みw2を掛け合わせてmatmul重み付けします。この値がモデルの予測値y_predとなります。

h = tf.matmul(x, w1)
h_relu = tf.maximum(h, tf.zeros(1))
y_pred = tf.matmul(h_relu, w2)

8. 損失の計算

モデルが予測した値y_predと答えyとの間の二乗誤差を計算しこれを損失lossとします。

tf.reduce_sumは、テンソルの足し算総和の計算をします。

loss = tf.reduce_sum((y - y_pred) *2.0)

9. 重み(Weight)の勾配計算

これより先は、モデルが予測した値y_predと答えyを見比べて、正しく答えyを予測できるようにモデルのパラメータを更新していきます。

具体的には、重みw1, w2の勾配(grad_w1, grad_w2)を計算します。

TensorFlowでは、gradientsメソッドで勾配が計算できます。

grad_w1, grad_w2 = tf.gradients(loss, [w1, w2])

10. 重み(w)の更新

計算した勾配(grad_w1, grad_w2)をもとに、重み(w1, w2)を更新します。自動微分を用いた場合のパラメータの更新は、torch.no_grad()で対象のパラメータをくくることで計算できます。

確率勾配降下法(SGD: stochastic gradient descent)は、重みを更新する上でよく使われる最適化アルゴリズムで、以下の式で表されます。

weight = weight - learning_rate gradient

今回は、学習率を1e-6、新しい重みw1, w2new_w1new_w2にして、SGDを用いて重みを更新します。

learning_rate = 1e-6
new_w1 = w1.assign(w1 - learning_rate grad_w1)
new_w2 = w2.assign(w2 - learning_rate grad_w2)

11. 計算グラフの実行

設計した計算グラフをこのコード部分で実行する。

with tf.Session() as sess:
    # Run the graph once to initialize the Variables w1 and w2.
    sess.run(tf.global_variables_initializer())

    # Create numpy arrays holding the actual data for the inputs x and targets
    # y
    x_value = np.random.randn(N, D_in)
    y_value = np.random.randn(N, D_out)
    for t in range(500):
        # Execute the graph many times.
        loss_value, _, _ = sess.run([loss, new_w1, new_w2],
                                    feed_dict={x: x_value, y: y_value})
        if t % 100 == 99:
            print(t, loss_value)

TensorFlowで、計算グラフを実行するには、with tf.Session() as sess:でSessionオブジェクトを作成する。

with tf.Session() as sess:

重みw1w2を初期化するため、Sessionをrunする。

    # Run the graph once to initialize the Variables w1 and w2.
    sess.run(tf.global_variables_initializer())

ニューラルネットワークに入力するデータx_valueとニューラルネットワークが予測すべき値y_valueを乱数np.random.randn`で決める。

    # Create numpy arrays holding the actual data for the inputs x and targets
    # y
    x_value = np.random.randn(N, D_in)
    y_value = np.random.randn(N, D_out)

今回は、学習回数を500回とする。

各学習ごとに、損失関数loss、更新された重みw1w2から損失値(二乗誤差)を計算する。

100回学習するごとに、学習回数tと二乗誤差loss_valueを出力する。

    for t in range(500):
        loss_value, _, _ = sess.run([loss, new_w1, new_w2],
                                    feed_dict={x: x_value, y: y_value})
        if t % 100 == 99:
            print(t, loss_value)

12. 実行

以下のコードを5_static_graphs.pyとして保存します。

12.1. 5_static_graphs.py

import tensorflow as tf
import numpy as np

# First we set up the computational graph:

# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

# Create placeholders for the input and target data; these will be filled
# with real data when we execute the graph.
x = tf.placeholder(tf.float32, shape=(None, D_in))
y = tf.placeholder(tf.float32, shape=(None, D_out))

# Create Variables for the weights and initialize them with random data.
# A TensorFlow Variable persists its value across executions of the graph.
w1 = tf.Variable(tf.random_normal((D_in, H)))
w2 = tf.Variable(tf.random_normal((H, D_out)))

# Forward pass
h = tf.matmul(x, w1)
h_relu = tf.maximum(h, tf.zeros(1))
y_pred = tf.matmul(h_relu, w2)

# Compute loss using operations on TensorFlow Tensors
loss = tf.reduce_sum((y - y_pred) *2.0)

# Compute gradient of the loss with respect to w1 and w2.
grad_w1, grad_w2 = tf.gradients(loss, [w1, w2])

# Update the weights using gradient descent.
learning_rate = 1e-6
new_w1 = w1.assign(w1 - learning_rate grad_w1)
new_w2 = w2.assign(w2 - learning_rate grad_w2)

# actually execute the graph.
with tf.Session() as sess:
    # Run the graph once to initialize the Variables w1 and w2.
    sess.run(tf.global_variables_initializer())

    # Create numpy arrays holding the actual data for the inputs x and targets
    # y
    x_value = np.random.randn(N, D_in)
    y_value = np.random.randn(N, D_out)
    for t in range(500):
        loss_value, _, _ = sess.run([loss, new_w1, new_w2],
                                    feed_dict={x: x_value, y: y_value})
        if t % 100 == 99:
            print(t, loss_value)

保存ができたら実行しましょう。

左の数字が学習回数、右の数値がパーモデルの推定値と実際の答えと二乗誤差です。

学習を重ねるごとに、二乗誤差が小さくなることがわかります。

$ python3 5_static_graphs.py
99 1597.0376
199 16.325699
299 0.24026018
399 0.0046504317
499 0.00029210464

【PyTorch】サンプル④ 〜Defining New autograd Functions(自動微分関数の定義)〜


1. 目的
2. 前準備
3. PyTorchのインポート
4. クラスの定義
5. データタイプとデバイスの選択
6. 使用するデータ
7. 重み付けの初期化
8. 学習
9. データの入力 (forward)
10. 損失の計算
11. 重み(Weight)の勾配計算
12. 重みの更新
13. 実行
13.1. 4_autograd_func.py


1. 目的

PyTorchのチュートリアルPyTorch: Defining New autograd Functionsを参考にPyTorchテンソル(tensor)と自動微分(autograd)を使って、損失(loss)や重み(weight)の計算をする。

前回の、【PyTorch】サンプル③ 〜TENSORS AND AUTOGRAD(テンソルと自動微分)〜では、自動微分の基本的な扱いをご紹介しました。

今回は、前回使用した自動微分のコードを関数化してみましょう。

2. 前準備

PyTorchのインストールはこちらから。

初めて、Google Colaboratoryを使いたい方は、こちらをご覧ください。

3. PyTorchのインポート

import torch

4. クラスの定義

自動微分のコード部分を以下のようなMyReLUというクラスで定義します。

class MyReLU(torch.autograd.Function):

    @staticmethod
    def forward(ctx, input):
        ctx.save_for_backward(input)
        return input.clamp(min=0)

    @staticmethod
    def backward(ctx, grad_output):
        input, = ctx.saved_tensors
        grad_input = grad_output.clone()
        grad_input[input < 0] = 0
        return grad_input

ここでは、torch.autograd.Functionを継承したサブクラスでMyReLUを定義します。

class MyReLU(torch.autograd.Function):

こちらが、パーセプトロンにデータをインプット(forward)する部分です。

@staticmethodはスタティックメソッドといって、そのクラスでしか使用しない関数を定義する時に使います。

また、スタティックメソッドの特徴として、クラスでよくあるselfを第一引数として取りません。

ctxは、逆伝播(backward)計算の際に用いる隠しオブジェクトです。

ctx.save_for_backwardメソッドを使うことで逆伝播計算に用いるオブジェクトのキャッシュを取って置くことができる。

input.clamp(min=0)で、入力されたinputの要素の中で0以下は0、0以上はそのままの要素の値として出力を返す。

    @staticmethod
    def forward(ctx, input):
        ctx.save_for_backward(input)
        return input.clamp(min=0)

この部分が、逆伝播計算の部分。

ctx.saved_tensorsで、パーセプトロンに入力された値inputを取り出します。

パーセプトロンからの出力勾配をgrad_output.clone()でコピーする。

    @staticmethod
    def backward(ctx, grad_output):
        input, = ctx.saved_tensors
        grad_input = grad_output.clone()
        grad_input[input < 0] = 0
        return grad_input
`

このように、PyTorchにおける自動微分関数の定義では、`forward`と`backward`をクラスのメソッドとして定義します。

###  5. <a name='-3'></a>データタイプとデバイスの選択
これから作成するテンソルのデータ型`dtype`を`float`にします。

テンソル計算をするデバイスは、`torch.device`で指定することができます。
今回は、CPUを使用することにします。
[python]
dtype = torch.float
device = torch.device("cpu")

GPUを使う方は、deviceの部分を以下のコードに変えて下さい。

device = torch.device("cuda:0") # Uncomment this to run on GPU

6. 使用するデータ

バッチサイズNを64、入力の次元D_inを1000、隠れ層の次元Hを100、出力の次元D_outを10とします。

# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

入力(x)と予測したい(y)を乱数で定義します。

PyTorchテンソルを定義する場合には、どのデバイスでdevice、どのデータ型dtypeかを指定することができます。

# Create random input and output data
x = torch.randn(N, D_in, device=device, dtype=dtype)
y = torch.randn(N, D_out, device=device, dtype=dtype)

7. 重み付けの初期化

乱数を使って重みを初期化します。

自動微分をしたい場合には、初期化する時点でrequires_gradのフラグをTrueにしておきます。デフォルトはrequires_grad=Falseです。

# Randomly initialize weights
w1 = torch.randn(D_in, H, device=device, dtype=dtype, requires_grad=True)
w2 = torch.randn(H, D_out, device=device, dtype=dtype, requires_grad=True)

8. 学習

学習率を1e-6として、学習回数を500回とします。

learning_rate = 1e-6
for t in range(500):

9. データの入力 (forward)

データをパーセプトロンに入力して、パーセプトロンの予測値y_predを計算します。

以前の自動微分のサンプルでは、

    y_pred = x.mm(w1).clamp(min=0).mm(w2)

このように記述しましたが、今回は自動微分の関数をMyReLUとしてに定義してあるので、上のコードをMyReLUを用いて以下のように書き換えます。

    # Forward pass: compute predicted y
    relu = MyReLU.apply
    y_pred = relu(x.mm(w1)).mm(w2)

10. 損失の計算

パーセプトロンが予測した値(y_pred)と答え(y)との間の二乗誤差を計算しこれを損失(loss)とします。

損失の値を、tensor配列からスカラー値で取得したい場合には、item()メソッドを用います。

各学習回数ごとに、学習回数(t)と二乗誤差(loss)を表示します。

今回は、結果を見やすくするためにif t % 100 == 99:の部分で、学習回数100回ごとに結果を出力します。
やっていることは、学習回数(t)を100で割ったときのあまりが99であるときにprint(t, loss)を実行です。

    # Compute and print loss
    loss = (y_pred - y).pow(2).sum()
    if t % 100 == 99:
        print(t, loss.item())

11. 重み(Weight)の勾配計算

これより先は、パーセプトロンが予測した値(y_pred)と答え(y)を見比べて、正しく答え(y)を予測できるようにパーセプトロンのパラメータを更新していきます。

    loss.backward()

12. 重みの更新

計算した勾配(grad_w1, grad_w2)をもとに、重み(w1, w2)を更新します。自動微分を用いた場合のパラメータの更新は、torch.no_grad()で対象のパラメータをくくることで計算できます。

確率勾配降下法(SGD: stochastic gradient descent)は、重みを更新する上でよく使われる最適化アルゴリズムで、以下の式で表されます。

weight = weight - learning_rate gradient

SGDを用いてパラメータを更新するには、このように記述します。

一度更新した、パラメータはgrad.zero()でゼロに初期化します。

    with torch.no_grad():
        w1 -= learning_rate w1.grad
        w2 -= learning_rate w2.grad

        # Manually zero the gradients after updating weights
        w1.grad.zero_()
        w2.grad.zero_()

13. 実行

以下のコードを4_autograd_func.pyとして保存します。

13.1. 4_autograd_func.py

import torch

dtype = torch.float
device = torch.device("cpu")
# device = torch.device("cuda:0") # Uncomment this to run on GPU

# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

# Create random Tensors to hold input and outputs.
x = torch.randn(N, D_in, device=device, dtype=dtype)
y = torch.randn(N, D_out, device=device, dtype=dtype)

# Create random Tensors for weights.
w1 = torch.randn(D_in, H, device=device, dtype=dtype, requires_grad=True)
w2 = torch.randn(H, D_out, device=device, dtype=dtype, requires_grad=True)

learning_rate = 1e-6
for t in range(500):
    # Forward pass: compute predicted y using operations on Tensors
    y_pred = x.mm(w1).clamp(min=0).mm(w2)

    # Compute and print loss using operations on Tensors.
    loss = (y_pred - y).pow(2).sum()
    if t % 100 == 99:
        print(t, loss.item())

    # Use autograd to compute the backward pass.
    loss.backward()

    # Manually update weights using gradient descent. Wrap in torch.no_grad()
    with torch.no_grad():
        w1 -= learning_rate w1.grad
        w2 -= learning_rate w2.grad

        # Manually zero the gradients after updating weights
        w1.grad.zero_()
        w2.grad.zero_()

保存ができたら実行しましょう。

左の数字が学習回数、右の数値がパーセプトロンの推定値と実際の答えと二乗誤差です。

学習を重ねるごとに、二乗誤差が小さくなることがわかります。

$ python3 4_autograd_func.py 
99 763.6893310546875
199 9.174347877502441
299 0.22351056337356567
399 0.007267479319125414
499 0.0004901751526631415

Ubuntu 20.04 の mini.iso

Ubuntu 20.04 の mini.iso はかなり隠れたところにあります。
こちらから入手できます。

http://archive.ubuntu.com/ubuntu/dists/focal/main/installer-amd64/current/legacy-images/netboot/

日本のミラーサイトはこちら。
http://jp.archive.ubuntu.com/ubuntu/dists/focal/main/installer-amd64/current/legacy-images/netboot/

いつも探すので備忘録として書いておきます。

【PyTorch】サンプル③ 〜TENSORS AND AUTOGRAD(テンソルと自動微分)〜


1. 目的
2. 前準備
3. 自動微分(AUTOGRAD)
4. PyTorchのインポート
5. データタイプとデバイスの選択
6. 使用するデータ
7. 重み付けの初期化
8. 学習
9. データの入力
10. 損失の計算
11. 重み(Weight)の勾配計算
12. 重みの更新
13. 実行
13.1. 3_autograd.py


1. 目的

PyTorchのチュートリアルPyTorch: Tensors and autogradを参考にPyTorchテンソル(tensor)と自動微分(autograd)を使って、損失(loss)や重み(weight)の計算をする。

これまでは、PyTorchに実装されている自動微分機能を使わずにニューラルネットワークのパラメータの勾配を計算していましたが、PyTorchの自動微分(autograd)機能を使うとパラメータの勾配計算簡単にすることができます。

2. 前準備

PyTorchのインストールはこちらから。

初めて、Google Colaboratoryを使いたい方は、こちらをご覧ください。

3. 自動微分(AUTOGRAD)

PyTorchにおけるテンソルは計算グラフでは、ノード(node)と表現されます。

例えば、テンソルxx.requires_grad=Trueとフラグを設定すると、x.gradでxの勾配を計算することができます。

4. PyTorchのインポート

import torch

5. データタイプとデバイスの選択

これから作成するテンソルのデータ型dtypefloatにします。

テンソル計算をするデバイスは、torch.deviceで指定することができます。
今回は、CPUを使用することにします。

dtype = torch.float
device = torch.device("cpu")

GPUを使う方は、deviceの部分を以下のコードに変えて下さい。

device = torch.device("cuda:0") # Uncomment this to run on GPU

6. 使用するデータ

バッチサイズNを64、入力の次元D_inを1000、隠れ層の次元Hを100、出力の次元D_outを10とします。

# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

入力(x)と予測したい(y)を乱数で定義します。

PyTorchテンソルを定義する場合には、どのデバイスでdevice、どのデータ型dtypeかを指定することができます。

# Create random input and output data
x = torch.randn(N, D_in, device=device, dtype=dtype)
y = torch.randn(N, D_out, device=device, dtype=dtype)

7. 重み付けの初期化

乱数を使って重みを初期化します。

注目してほしい点は、requires_grad=Trueです。
自動微分をしたい場合には、初期化する時点でrequires_gradのフラグをTrueにしておきます。デフォルトはrequires_grad=Falseです。

# Randomly initialize weights
w1 = torch.randn(D_in, H, device=device, dtype=dtype, requires_grad=True)
w2 = torch.randn(H, D_out, device=device, dtype=dtype, requires_grad=True)

8. 学習

学習率を1e-6として、学習回数を500回とします。

learning_rate = 1e-6
for t in range(500):

9. データの入力

二次元テンソルの入力(x)と重み(w1)をmmで掛け算することで重み付けをします(h)。(一次元テンソルの積は、dotで計算できます。)

重み付けした値の要素から、clamp(min=~0)で、0以上の要素は残し、0以下のものは0とします。

最後に、重み(w2)を掛け合わせて重み付けしますmm(w2)
この値がパーセプトロンの予測値(y_pred)となります。

    # Forward pass: compute predicted y
    y_pred = x.mm(w1).clamp(min=0).mm(w2)

10. 損失の計算

パーセプトロンが予測した値(y_pred)と答え(y)との間の二乗誤差を計算しこれを損失(loss)とします。

損失の値を、tensor配列からスカラー値で取得したい場合には、item()メソッドを用います。

各学習回数ごとに、学習回数(t)と二乗誤差(loss)を表示します。

今回は、結果を見やすくするためにif t % 100 == 99:の部分で、学習回数100回ごとに結果を出力します。
やっていることは、学習回数(t)を100で割ったときのあまりが99であるときにprint(t, loss)を実行です。

    # Compute and print loss
    loss = (y_pred - y).pow(2).sum()
    if t % 100 == 99:
        print(t, loss.item())

11. 重み(Weight)の勾配計算

これより先は、パーセプトロンが予測した値(y_pred)と答え(y)を見比べて、正しく答え(y)を予測できるようにパーセプトロンのパラメータを更新していきます。

NumPyを用いたサンプルPyTorchを用いたサンプルでは、ガリガリ勾配の計算式を書いていましたが、自動微分を使えば、一行で勾配を計算できます

    loss.backward()

12. 重みの更新

計算した勾配(grad_w1, grad_w2)をもとに、重み(w1, w2)を更新します。自動微分を用いた場合のパラメータの更新は、torch.no_grad()で対象のパラメータをくくることで計算できます。

確率勾配降下法(SGD: stochastic gradient descent)は、重みを更新する上でよく使われる最適化アルゴリズムで、以下の式で表されます。

weight = weight - learning_rate gradient

SGDを用いてパラメータを更新するには、このように記述します。

一度更新した、パラメータはgrad.zero()でゼロに初期化します。

    with torch.no_grad():
        w1 -= learning_rate w1.grad
        w2 -= learning_rate w2.grad

        # Manually zero the gradients after updating weights
        w1.grad.zero_()
        w2.grad.zero_()

13. 実行

以下のコードを3_autograd.pyとして保存します。

13.1. 3_autograd.py

import torch

dtype = torch.float
device = torch.device("cpu")
# device = torch.device("cuda:0") # Uncomment this to run on GPU

# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

# Create random Tensors to hold input and outputs.
x = torch.randn(N, D_in, device=device, dtype=dtype)
y = torch.randn(N, D_out, device=device, dtype=dtype)

# Create random Tensors for weights.
w1 = torch.randn(D_in, H, device=device, dtype=dtype, requires_grad=True)
w2 = torch.randn(H, D_out, device=device, dtype=dtype, requires_grad=True)

learning_rate = 1e-6
for t in range(500):
    # Forward pass: compute predicted y using operations on Tensors
    y_pred = x.mm(w1).clamp(min=0).mm(w2)

    # Compute and print loss using operations on Tensors.
    loss = (y_pred - y).pow(2).sum()
    if t % 100 == 99:
        print(t, loss.item())

    # Use autograd to compute the backward pass.
    loss.backward()

    # Manually update weights using gradient descent. Wrap in torch.no_grad()
    with torch.no_grad():
        w1 -= learning_rate w1.grad
        w2 -= learning_rate w2.grad

        # Manually zero the gradients after updating weights
        w1.grad.zero_()
        w2.grad.zero_()

保存ができたら実行しましょう。

左の数字が学習回数、右の数値がパーセプトロンの推定値と実際の答えと二乗誤差です。

学習を重ねるごとに、二乗誤差が小さくなることがわかります。

$ python3 3_autograd.py 
99 743.60400390625
199 7.154838562011719
299 0.13969427347183228
399 0.0036201069597154856
499 0.000234781633480452

【PyTorch】サンプル② 〜TENSOR (テンソル) 〜


1. 目的
2. 前準備
3. TENSOR(テンソル)
4. PyTorchのインポート
5. データタイプとデバイスの選択
6. 使用するデータ
7. 重み付けの初期化
8. 学習
9. データの入力
10. 損失の計算
11. 重み(Weight)の勾配計算
12. 重みの更新
13. 実行
13.1. 2_tensor.py


### 1. 目的
PyTorchのチュートリアルPyTorch: Tensorsを参考にPyTorch Tensor(テンソル)を使って、損失(loss)や重み(weight)の計算をする。

PyTorchの特徴の一つである、テンソルとNumPyの違いを手を動かしながら感じる。

NumPyを用いた例は、こちらから。

2. 前準備

PyTorchのインストールはこちらから。

初めて、Google Colaboratoryを使いたい方は、こちらをご覧ください。

3. TENSOR(テンソル)

PyTorchのtensor(テンソル)は、基本的にはnumpy配列と同じです。

しかし、PyTorchのテンソルは、deep learning(深層学習)や計算グラフ、勾配といった概念に合わせた配列です。

PyTorchテンソルだけがCPUとGPUのどちらでも計算が実行でき、これが、numpy配列とPyTorchテンソルの最も大きな違いです。

4. PyTorchのインポート

import torch

5. データタイプとデバイスの選択

これから作成するテンソルのデータ型dtypefloatにします。

テンソル計算をするデバイスは、torch.deviceで指定することができます。
今回は、CPUを使用することにします。

dtype = torch.float
device = torch.device("cpu")

GPUを使う方は、deviceの部分を以下のコードに変えて下さい。

device = torch.device("cuda:0") # Uncomment this to run on GPU

6. 使用するデータ

バッチサイズNを64、入力の次元D_inを1000、隠れ層の次元Hを100、出力の次元D_outを10とします。

# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

入力(x)と予測したい(y)を乱数で定義します。

PyTorchテンソルを定義する場合には、どのデバイスでdevice、どのデータ型dtypeかを指定することができます。

# Create random input and output data
x = torch.randn(N, D_in, device=device, dtype=dtype)
y = torch.randn(N, D_out, device=device, dtype=dtype)

7. 重み付けの初期化

乱数を使って重みを初期化します。

# Randomly initialize weights
w1 = torch.randn(D_in, H, device=device, dtype=dtype)
w2 = torch.randn(H, D_out, device=device, dtype=dtype)

8. 学習

学習率を1e-6として、学習回数を500回とします。

learning_rate = 1e-6
for t in range(500):

9. データの入力

二次元テンソルの入力(x)と重み(w1)をmmで掛け算することで重み付けをします(h)。(一次元テンソルの積は、dotで計算できます。)

重み付けした値(h)の要素から、clamp(min=~0)で、0以上のものは残し、0以下のものは0とします。numpyでは、numpy.maximum(h,0)と書きます。

最後に、重み(w2)を掛け合わせて重み付けしますmm(w2)
この値がパーセプトロンの予測値(y_pred)となります。

    # Forward pass: compute predicted y
    h = x.mm(w1)
    h_relu = h.clamp(min=0)
    y_pred = h_relu.mm(w2)

10. 損失の計算

パーセプトロンが予測した値(y_pred)と答え(y)との間の二乗誤差を計算しこれを損失(loss)とします。
numpyでは、二乗計算をnumpy.squareで実行していましたが、PyTorchテンソルでの二乗計算はpow(2)と記述します。

損失の値を、tensor配列からスカラー値で取得したい場合には、item()メソッドを用います。

各学習回数ごとに、学習回数(t)と二乗誤差(loss)を表示します。

今回は、結果を見やすくするためにif t % 100 == 99:の部分で、学習回数100回ごとに結果を出力します。
やっていることは、学習回数(t)を100で割ったときのあまりが99であるときにprint(t, loss)を実行です。

    # Compute and print loss
    loss = (y_pred - y).pow(2).sum().item()
    if t % 100 == 99:
        print(t, loss)

11. 重み(Weight)の勾配計算

これより先は、パーセプトロンが予測した値(y_pred)と答え(y)を見比べて、正しく答え(y)を予測できるようにパーセプトロンのパラメータを更新していきます。

まず、重み(w1, w2)の勾配(grad_w1, grad_w2)を計算します。

numpyと違う点は、表のとおりです。

PyTorch NumPy
配列の転置 .t() .T
二次元配列の掛け算 .mm() .dot()
配列のコピー .clone() .copy()
    grad_y_pred = 2.0 (y_pred - y)
    grad_w2 = h_relu.t().mm(grad_y_pred)
    grad_h_relu = grad_y_pred.mm(w2.t())
    grad_h = grad_h_relu.clone()
    grad_h[h < 0] = 0
    grad_w1 = x.t().mm(grad_h)

12. 重みの更新

計算した勾配(grad_w1, grad_w2)をもとに、重み(w1, w2)を更新します。

確率勾配降下法(SGD: stochastic gradient descent)は、重みを更新する上でよく使われる最適化アルゴリズムで、以下の式で表されます。
weight = weight - learning_rate gradient

SGDは、以下のコードで実行できます。
ここの部分は、numpyもPyTorchテンソルも同じコードです。

    # Update weights using gradient descent
    w1 -= learning_rate grad_w1
    w2 -= learning_rate grad_w2

13. 実行

以下のコードを2_tensor.pyとして保存します。

13.1. 2_tensor.py

import torch


dtype = torch.float
device = torch.device("cpu")
# device = torch.device("cuda:0") # Uncomment this to run on GPU

# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

# Create random input and output data
x = torch.randn(N, D_in, device=device, dtype=dtype)
y = torch.randn(N, D_out, device=device, dtype=dtype)

# Randomly initialize weights
w1 = torch.randn(D_in, H, device=device, dtype=dtype)
w2 = torch.randn(H, D_out, device=device, dtype=dtype)

learning_rate = 1e-6
for t in range(500):
    # Forward pass: compute predicted y
    h = x.mm(w1)
    h_relu = h.clamp(min=0)
    y_pred = h_relu.mm(w2)

    # Compute and print loss
    loss = (y_pred - y).pow(2).sum().item()
    if t % 100 == 99:
        print(t, loss)

    # Backprop to compute gradients of w1 and w2 with respect to loss
    grad_y_pred = 2.0 (y_pred - y)
    grad_w2 = h_relu.t().mm(grad_y_pred)
    grad_h_relu = grad_y_pred.mm(w2.t())
    grad_h = grad_h_relu.clone()
    grad_h[h < 0] = 0
    grad_w1 = x.t().mm(grad_h)

    # Update weights using gradient descent
    w1 -= learning_rate grad_w1
    w2 -= learning_rate grad_w2

保存ができたら実行しましょう。

左の数字が学習回数、右の数値がパーセプトロンの推定値と実際の答えと二乗誤差です。

学習を重ねるごとに、二乗誤差が小さくなることがわかります。

$ python3 2_tensor.py 
99 306.79400634765625
199 0.4559670686721802
299 0.0014809096464887261
399 6.634114834014326e-05
499 1.7021899111568928e-05

【PyTorch】サンプル① 〜NUMPY〜


1. 目的
2. 前準備
3. 予備知識
4. NumPyのインポート
5. データ
6. 重み付けの初期化
7. 学習
8. データの入力
9. 重み(Weight)の勾配計算
10. 損失の計算
11. 重みの更新
12. 実行
12.1. 1_numpy.py


1. 目的

PyTorchのチュートリアルWarm-up: numpyを参考にNumpyを使って、損失(loss)や重み(weight)の計算をする。

PyTorchの特徴の一つである、テンソルとNumpyの違いを理解するための前準備。

2. 前準備

PyTorchのインストールはこちらから。

初めて、Google Colaboratoryを使いたい方は、こちらをご覧ください。

3. 予備知識

脳神経細胞は、樹上突起(Dendrites)、細胞体(Soma)、核(Nucleus)、軸索(Axon)、軸索終末(Axon Terminals)で構成されます。

(Credit: https://commons.wikimedia.org/wiki/File:Neuron_-_annotated.svg)

この脳神経細胞を数学的にモデルしたのが、パーセプトロンです。
神経樹上から細胞体に向けてやってくる信号(Xn)は、すべてが重要であるわけではありません。
入力される信号の中には、必要なものとそうでないものが混じっていると考えて、各信号に対して重み付け(Wn)をします。
神経樹上からの信号(Xn)は重みづけ(Wn)され、合算(Sum Σ)されます。
その合算した信号(z)が、その神経細胞を活性化させるかどうかを活性化関数(Activation function σ)でモデルします。
最後に、活性化関数からの出力に対して重み付けをして次のパーセプトロンに信号(a)を渡します。

(Credit: https://pythonmachinelearning.pro/perceptrons-the-first-neural-networks/)

今回は、パーセプトロンを用いて、入力される信号(x)からyを予測するケースを考えます。

4. NumPyのインポート

import numpy as np

5. データ

バッチサイズNを64、入力の次元D_inを1000、隠れ層の次元Hを100、出力の次元D_outを10とします。

# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

入力(x)と予測したい(y)を乱数で定義します。

# Create random input and output data
x = np.random.randn(N, D_in)
y = np.random.randn(N, D_out)

6. 重み付けの初期化

乱数を使って重みを初期化します。

# Randomly initialize weights
w1 = np.random.randn(D_in, H)
w2 = np.random.randn(H, D_out)

7. 学習

学習率を1e-6として、学習回数を500回とします。

learning_rate = 1e-6
for t in range(500):

8. データの入力

入力(x)と重み(w1)を掛け算.dotすることで重み付けをします(h)。
重み付けした値(h)の要素から、np.maximum(h,0)で、0以上のものは残し、0以下のものは0とします。
最後に、重み(w2)を掛け合わせて重み付けします。この値がパーセプトロンの予測値(y_pred)となります。

    # Forward pass: compute predicted y
    h = x.dot(w1)
    h_relu = np.maximum(h, 0)
    y_pred = h_relu.dot(w2)

9. 重み(Weight)の勾配計算

これより先は、パーセプトロンが予測した値(y_pred)と答え(y)を見比べて、正しく答え(y)を予測できるようにパーセプトロンのパラメータを更新していきます。

まず、重み(w1, w2)の勾配(grad_w1, grad_w2)を計算します。

    # Backprop to compute gradients of w1 and w2 with respect to loss
    grad_y_pred = 2.0 * (y_pred - y)
    grad_w2 = h_relu.T.dot(grad_y_pred)
    grad_h_relu = grad_y_pred.dot(w2.T)
    grad_h = grad_h_relu.copy()
    grad_h[h < 0] = 0
    grad_w1 = x.T.dot(grad_h)

10. 損失の計算

パーセプトロンが予測した値(y_pred)と答え(y)との間の二乗誤差を計算しこれを損失(loss)とします。
np.squreareでy_predとyの差を二乗して、sum()で平均しています。
各学習回数ごとに、学習回数(t)と二乗誤差(loss)を表示します。

    # Compute and print loss
    loss = np.square(y_pred - y).sum()
    print(t, loss)

11. 重みの更新

計算した勾配(grad_w1, grad_w2)をもとに、重み(w1, w2)を更新します。

確率勾配降下法(SGD: stochastic gradient descent)は、重みを更新する上でよく使われる最適化アルゴリズムで、以下の式で表されます。

weight = weight - learning_rate * gradient

SGDは、以下のコードで実行できます。

    # Update weights
    w1 -= learning_rate * grad_w1
    w2 -= learning_rate * grad_w2

12. 実行

以下のコードを1_numpy.pyとして保存します。

12.1. 1_numpy.py

import numpy as np

# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

# Create random input and output data
x = np.random.randn(N, D_in)
y = np.random.randn(N, D_out)

# Randomly initialize weights
w1 = np.random.randn(D_in, H)
w2 = np.random.randn(H, D_out)

learning_rate = 1e-6
for t in range(500):
    # Forward pass: compute predicted y
    h = x.dot(w1)
    h_relu = np.maximum(h, 0)
    y_pred = h_relu.dot(w2)

    # Compute and print loss
    loss = np.square(y_pred - y).sum()
    print(t, loss)

    # Backprop to compute gradients of w1 and w2 with respect to loss
    grad_y_pred = 2.0 * (y_pred - y)
    grad_w2 = h_relu.T.dot(grad_y_pred)
    grad_h_relu = grad_y_pred.dot(w2.T)
    grad_h = grad_h_relu.copy()
    grad_h[h < 0] = 0
    grad_w1 = x.T.dot(grad_h)

    # Update weights
    w1 -= learning_rate * grad_w1
    w2 -= learning_rate * grad_w2

保存ができたら実行しましょう。
左の数字が学習回数、右の数値がパーセプトロンの推定値と実際の答えと二乗誤差です。
学習を重ねるごとに、二乗誤差が小さくなることがわかります。

$ python3 1_numpy.py 
0 33318410.89325847
1 33449484.266180404
2 42189212.89431849
3 51379306.420906566
4 48992878.8013583

...

499 1.529654790074364e-05