BashコマンドとPythonの対応表

BashコマンドとPythonの対応表

私は bash などのシェルでのディレクトリ操作、ファイル操作の頻用コマンドは頭に入っているのですが、Pythonで同じことをやろうとするといつも戸惑います。なので、一覧表を作成してみました。なお、近年は os モジュールよりも pathlib モジュールがよく使われるとのことでその対応表もまとめてみました。なお、 os モジュールの説明中に shutil も入っています。コピーや削除はこちらを使うので、ここに入れておきます。

続きを読む

【Python】DICOMヘッダーをCSVに保存

1. 目的

PythonのPydicomライブラリを用いて、DICOMヘッダーをCSVにまとめて保存

2. 準備

2.1. ライブラリの準備

Pydicomは、DICOMのヘッダーや画像を操作するのに用いるライブラリである。

Pydicomのインストールは、以下のコマンドを実行。

pip3 install pydicom

CSV形式の表データを扱うには、Pandasライブラリを用いる。

Pandasライブラリのインストールは、以下のコマンドを実行。

pip3 install pandas

2.2. データの準備

次のような、フォルダ構造でデータを準備する。この場合では、各被験者フォルダの中にDICOMが保存されている。

DICOM_folder
├── Subject001
│   ├── XXX.dcm
│   ├── ...
│   └── XXX.dcm
├── Subject002
│   ├── XXX.dcm
│   ├── ...
│   └── XXX.dcm
├── ...
└── SubjectXXX

2.3. スクリプトの準備

次のコードを、extract_dcm_header.pyとして保存する。このとき、スクリプトはDICOM_folderフォルダと同じ階層に保存する。

import os
import pydicom
import pandas as pd

input='DICOM_folder'  # Input folder
output='dicom_headers.csv'  # Output CSV

dcm_dfs = []
failed_files = []
processed_files = []
for root, _, files in os.walk(input):  # Find DICOM file for each subject
    if len(files) != 0:  # If DICOM files exist
        try:
            f = os.path.join(root, files[0])
            dcm = pydicom.dcmread(f)  # Read DICOM
            _df = pd.DataFrame({dcm[k].keyword: [dcm[k].value] for k in dcm.keys() if dcm[k].keyword != "PixelData"})  # Read Headers
            dcm_dfs.append(_df)  # Gather headers of all subjects in a list
            processed_files.append(f)
        except:
            failed_files.append(f)

dcm_dfs = pd.concat(dcm_dfs, ignore_index=True)  # Concat headers of all subjects in a table
dcm_dfs.to_csv(output, index=False)  # Save as CSV

3. プログラムの実行

2.3. スクリプトの準備で用意した、extract_dcm_header.pyを実行するには、次のコマンドを実行する。

python3 ./extract_dcm_header.py

4. 結果の確認

収集したDICOMヘッダーは、dicom_headers.csvとして保存される。

5. コードの解説

まず、必要なライブラリを読み込む。

import os
import pydicom
import pandas as pd

ここでは、入力となるDICOMフォルダーと出力となるDICOMヘッダーのまとまったCSVの名前を定義している。

今回の場合だとinput='DICOM_folder'output='dicom_headers.csv'

input='DICOM_folder'  # Input folder
output='dicom_headers.csv'  # Output CSV

データを格納するための、箱(リスト)を定義。

dcm_dfs = []
failed_files = []
processed_files = []

被験者ごとのDICOMファイルを検索。

for root, _, files in os.walk(input):  # Find DICOM file for each subject

DICOMファイルがある場合のみ、処理を実行。

    if len(files) != 0:  # If DICOM files exist

Pydicomを用いて、DICOMデータを読み込む。

        try:
            f = os.path.join(root, files[0])
            dcm = pydicom.dcmread(f)  # Read DICOM

DICOMからヘッダー(Header)情報を、Pandasで読み込む。

PixelDataタグを含めると、出力(CSV)が崩れておかしくなるので、収集に含めないようにしている。

            _df = pd.DataFrame({dcm[k].keyword: [dcm[k].value] for k in dcm.keys() if dcm[k].keyword != "PixelData"})  # Read Headers

収集した結果を、被験者ごとに処理をして、一つの箱(リスト)にまとめる。

            dcm_dfs.append(_df)  # Gather headers of all subjects in a list
            processed_files.append(f)
        except:
            failed_files.append(f)

すべての被験者のヘッダー情報を、一つの表形式のデータ(DataFrame型)に変換する。


dcm_dfs = pd.concat(dcm_dfs, ignore_index=True)  # Concat headers of all subjects in a table

結果を、CSVとして保存する。

dcm_dfs.to_csv(output, index=False)  # Save as CSV

Pythonでオブジェクトの型に準備されている通常メソッドの一覧を出力する関数

Pythonを勉強していると、「この型のメソッドは何だろう?」と思う時があります。
この時、オブジェクトを obj とすると

dir(obj)

とすることで、一覧を得ることができます。

たとえば、リスト型のメソッドを知りたいとします。dir() を使うと以下のようになります。

x = [1, 2]
dir(x)

[code lang=text]
['__add__', '__class__', '__class_getitem__',
 '__contains__', '__delattr__', '__delitem__',
 '__dir__', '__doc__', '__eq__', '__format__',
 '__ge__', '__getattribute__', '__getitem__',
 '__gt__', '__hash__', '__iadd__', '__imul__',
 '__init__', '__init_subclass__', '__iter__',
 '__le__', '__len__', '__lt__', '__mul__',
 '__ne__', '__new__', '__reduce__', '__reduce_ex__',
 '__repr__', '__reversed__', '__rmul__', '__setattr__',
 '__setitem__', '__sizeof__', '__str__', '__subclasshook__',
 'append', 'clear', 'copy', 'count', 'extend',
 'index', 'insert', 'pop', 'remove', 'reverse', 'sort']
[/code]

ここで、”–” からはじまるメソッドは特殊メソッドと言われ、その型の振る舞いを細かく調整するものとのことです。今回はここには踏み込みません。

今、私は、「通常メソッドだけリストアップしたい」と思いました。どうしたらできるでしょうか?

続きを読む

LinuxやmacOSでPythonの仮想環境を構築する方法

過去に、Anacondaに頼らない、pipとvenvを用いたPython環境の構築 という記事を書きました。今回、改めて、Pythonの仮想環境について理解が深まったので書きたいと思います。

仮想環境を構築したい背景

  • Python は、ひとつのシステムに様々なバージョンが存在しえます。macOSでの場合を、macOS の Python事情を理解するに解説しました。LinuxやWindowsも同じで複数のバージョンが存在しえます。
  • また、LinuxやmacOSにおいて、Pythonは、システムの重要なところを担っていたりします。Ubuntuであれば、 dpkg -l | grep python3 とすると、どれだけ多くの Python3に関連したパッケージがシステムにインストールされているかを確認することができます。

  • このような状況において、システムのPythonに追加でパッケージを入れていって、もし不具合が起きた場合、システムそのものが不安定になる可能性があります。

  • Pythonの仮想環境を使うと、システムの中に、独立したPythonの環境を構築することができます。「独立している」というのは、システムに一切影響を与えないということを意味します。不要になったらばっさり削除しても一切問題ありません。

  • そこで、以下で、仮想環境の構築の仕方を解説します

続きを読む

過去のPhilipsのfMRIデータをdcm2niixで変換しようとした時に4次元データにならない問題の解決法

ある施設のrs-fMRIのDICOMデータをNiftiに変換しようとした時に、以下のようになってしまい、4次元データができませんでした。

sub1_+rsfMRI_201.nii
sub1_+rsfMRI_201_t10000.nii
sub1_+rsfMRI_201_t100000.nii
sub1_+rsfMRI_201_t102500.nii
sub1_+rsfMRI_201_t105000.nii
sub1_+rsfMRI_201_t107500.nii
sub1_+rsfMRI_201_t110000.nii
sub1_+rsfMRI_201_t112500.nii
sub1_+rsfMRI_201_t115000.nii
sub1_+rsfMRI_201_t117500.nii
sub1_+rsfMRI_201_t120000.nii
sub1_+rsfMRI_201_t122500.nii
sub1_+rsfMRI_201_t12500.nii
...

ポイントは、ファイル名の後ろに tの後に数字がつくことです。

続きを読む

【AMICO】AMICOを用いた神経突起イメージング: NODDI



1. 目的
2. 準備
2.1. インストール
2.2. 使用データ
2.3. 前処理
3. 神経突起イメージング(NODDI)
3.1. AMICOのセットアップ
3.2. データの読み込み
3.3. 応答関数(response function)の計算
3.4. モデルフィッティング
3.5. 結果


1. 目的

  • AMICOを用いた神経突起イメージング: NODDI

2. 準備

2.1. インストール

Pythonを使って、AMICOを用いた神経突起イメージング(NODDI)をするために、以下のPythonパッケージをインストールする。

pip3 install dmri-amico

2.2. 使用データ

データを次のフォルダ構造で用意する。

Study/
└── Subject
    ├── DWI.nii.gz  # 拡散MRI
    ├── DWI_mask.nii.gz  # 拡散MRIマスク画像
    ├── bvals  # b-values
    └── bvecs  # b-vectors

2.3. 前処理

NODDI前に、拡散MRIの前処理をする。

  • 拡散MRIのノイズ除去(Software: MRtrix, DIPY)
  • ギブズのリンギングアーチファクト(Gibbs ringing)の除去(Software: MRtrix, DIPY)
  • 拡散MRIのバイアス(信号ムラ)補正(Software: MRtrix)
  • 拡散MRIの前処理 ~歪み・頭の動き・渦電流の補正(Software: FSL, MRtrix)

3. 神経突起イメージング(NODDI)

Pythonで以下のコマンドを実行。

3.1. AMICOのセットアップ

今回使用するファイル等の変数設定をする。

STUDY_DIR='Study'
SUBJECT_DIR='Subject'
DWI_FILE = 'DWI.nii.gz'
DWIMASK_FILE = 'DWI_mask.nii.gz'
BVALS_FILE = 'bvals'
BVECS_FILE = 'bvecs'

次に、使用するamicoパッケージのをインポートし、セットアップと初期化をする。

import amico
amico.core.setup()

3.2. データの読み込み

AMICOを実行するために、Study/Subjectフォルダ(PATH)を指定する。

ae = amico.Evaluation(STUDY_DIR, SUBJECT_DIR)

MPG軸情報(bvals/bvecs)の情報が入ったschemeファイルを生成する。

amico.util.fsl2scheme("{}/{}/{}".format(STUDY_DIR,SUBJECT_DIR,BVALS_FILE), "{}/{}/{}".format(STUDY_DIR,SUBJECT_DIR,BVECS_FILE),schemeFilename = "{}/{}/NODDI_protocol.scheme".format(STUDY_DIR,SUBJECT_DIR))
-> Writing scheme file to [ Study/Subject/NODDI_protocol.scheme ]
'Study/Subject/NODDI_protocol.scheme'

画像を読み込む。

ae.load_data(dwi_filename = DWI_FILE, scheme_filename = 'NODDI_protocol.scheme', mask_filename = DWIMASK_FILE, b0_thr = 0)
-> Loading data:
    * DWI signal
        - dim    = 130 x 130 x 82 x 129
        - pixdim = 1.769 x 1.769 x 1.800
    * Acquisition scheme
        - 129 samples, 2 shells
        - 1 @ b=0 , 64 @ b=1000.0 , 64 @ b=2000.0 
    * Binary mask
        - dim    = 130 x 130 x 82
        - pixdim = 1.769 x 1.769 x 1.800
        - voxels = 282878
   [ 4.4 seconds ]

-> Preprocessing:
    * Normalizing to b0... [ min=-3.28,  mean=0.25, max=22.86 ]
    * Keeping all b0 volume(s)
   [ 1.1 seconds ]

3.3. 応答関数(response function)の計算

NODDIモデルを設定して、応答関数(response function)を計算する。計算が完了するとkernelファイルが生成される。

ae.set_model('NODDI')
ae.generate_kernels()
-> Creating LUT for "NODDI" model:
   [==================================================] 100.0% 
   [ 373.3 seconds ]

作成したkernelファイルを読み込む。

ae.load_kernels()
-> Resampling LUT for subject "Subject":
   [==================================================] 100.0% 
   [ 112.8 seconds ]

3.4. モデルフィッティング

NODDIのモデルフィッティングを開始する。

ae.fit()
-> Fitting "NODDI" model to 282878 voxels:
   [==================================================] 100.0% 
   [ 02h 52m 07s ]

最後に、結果をNIfTIフォーマットで保存する。

ae.save_results()
-> Saving output to "AMICO/NODDI/*":
    - configuration  [OK]
    - FIT_dir.nii.gz  [OK]
    - FIT_ICVF.nii.gz  [OK]
    - FIT_OD.nii.gz  [OK]
    - FIT_ISOVF.nii.gz  [OK]
   [ DONE ]

3.5. 結果

結果は、「Study/Subject/AMICO/NODDI/」フォルダにある。

Study/Subject/AMICO/NODDI/
├── FIT_ICVF.nii.gz
├── FIT_ISOVF.nii.gz
├── FIT_OD.nii.gz
├── FIT_dir.nii.gz
└── config.pickle

画像はこちら。

【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

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


1. 目的
2. 準備
2.1. DIPYのインストール
2.2. 使用データ
3. 拡散MRIのノイズ除去
3.1. 必要なパッケージをインポート
3.2. 画像およびMPG軸情報の読み込み
3.3. マスク画像の作成
3.4. ギブズのリンギングアーチファクト除去
3.5. NIfTI形式で保存
3.6. 結果


### 1. 目的

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

2. 準備

2.1. DIPYのインストール

pip3 install dipy

2.2. 使用データ

データを次のフォルダ構造で用意する。

Study/
└── Subject
    ├── DWI.nii.gz  # 拡散MRI
    ├── DWI_mask.nii.gz  # 拡散MRIマスク画像
    ├── bvals  # b-values
    └── bvecs  # b-vectors

3. 拡散MRIのノイズ除去

Pythonで以下のコマンドを実行。

3.1. 必要なパッケージをインポート

from dipy.denoise.gibbs import gibbs_removal
import matplotlib.pyplot as plt
import numpy as np
from dipy.segment.mask import median_otsu
from dipy.io.image import load_nifti, save_nifti
from dipy.io.gradients import read_bvals_bvecs
from dipy.core.gradients import gradient_table

3.2. 画像およびMPG軸情報の読み込み

DWI_FILE = 'DWI.nii.gz'
BVALS_FILE = 'bvals'
BVECS_FILE = 'bvecs'

# Import data
data, affine = load_nifti(DWI_FILE)
bvals, bvecs = read_bvals_bvecs(BVALS_FILE, BVECS_FILE)
gtab = gradient_table(bvals, bvecs)

3.3. マスク画像の作成

median_otsu関数を用いて、b=0画像からマスク画像を生成する。vol_idxには、b0 volumeのvolume indexを渡す。

maskdata, mask = median_otsu(
    data, vol_idx=np.where(bvals == 0)[0])

3.4. ギブズのリンギングアーチファクト除去

gibbs_removal関数を用いて、リンギングアーチファクトを除去する。

data_corrected = gibbs_removal(maskdata)

3.5. NIfTI形式で保存

save_nifti関数で、画像をNIfTI形式で保存する。

save_nifti('DWI_degibbs.nii.gz', data_corrected.astype(np.float32), affine)

3.6. 結果

補正前後の画像は、以下の通り。

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

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

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

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

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

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

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

概要は以下になります。

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

続きを読む