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

PythonのPydicomライブラリを用いて、DICOMヘッダーをCSVにまとめて保存
Pydicomは、DICOMのヘッダーや画像を操作するのに用いるライブラリである。
Pydicomのインストールは、以下のコマンドを実行。
pip3 install pydicom
CSV形式の表データを扱うには、Pandasライブラリを用いる。
Pandasライブラリのインストールは、以下のコマンドを実行。
pip3 install pandas
次のような、フォルダ構造でデータを準備する。この場合では、各被験者フォルダの中にDICOMが保存されている。
DICOM_folder ├── Subject001 │ ├── XXX.dcm │ ├── ... │ └── XXX.dcm ├── Subject002 │ ├── XXX.dcm │ ├── ... │ └── XXX.dcm ├── ... └── SubjectXXX
次のコードを、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
2.3. スクリプトの準備で用意した、extract_dcm_header.pyを実行するには、次のコマンドを実行する。
python3 ./extract_dcm_header.py
収集したDICOMヘッダーは、dicom_headers.csvとして保存される。
まず、必要なライブラリを読み込む。
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を勉強していると、「この型のメソッドは何だろう?」と思う時があります。
この時、オブジェクトを 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]
ここで、”–” からはじまるメソッドは特殊メソッドと言われ、その型の振る舞いを細かく調整するものとのことです。今回はここには踏み込みません。
今、私は、「通常メソッドだけリストアップしたい」と思いました。どうしたらできるでしょうか?
過去に、Anacondaに頼らない、pipとvenvを用いたPython環境の構築 という記事を書きました。今回、改めて、Pythonの仮想環境について理解が深まったので書きたいと思います。
また、LinuxやmacOSにおいて、Pythonは、システムの重要なところを担っていたりします。Ubuntuであれば、 dpkg -l | grep python3 とすると、どれだけ多くの Python3に関連したパッケージがシステムにインストールされているかを確認することができます。
このような状況において、システムのPythonに追加でパッケージを入れていって、もし不具合が起きた場合、システムそのものが不安定になる可能性があります。
Pythonの仮想環境を使うと、システムの中に、独立したPythonの環境を構築することができます。「独立している」というのは、システムに一切影響を与えないということを意味します。不要になったらばっさり削除しても一切問題ありません。
そこで、以下で、仮想環境の構築の仕方を解説します
macOSでPythonを使おうとする時、様々な選択肢があります。
まず、それぞれのインストール方法とそのPythonのパスを明確にします。バージョンは2024年1月現在のものになります。
ある施設の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の後に数字がつくことです。
1. 目的
2. 準備
2.1. data.csv
3. ソースコード
4. 結果
4.1. パラメトリック検定の結果(para_result.csv)
4.2. パラメトリック検定の結果(nonpara_result.csv)

1. 目的
2. 準備
3. 必要なpackageのdownload
4. 必要なmoduleのimport
5. CSVファイルの読み込み
6. 2変数の定義
7. MAEの計算
8. RMSEの計算
9. 使用したコードまとめ

1. 目的
2. 準備
3. フォルダ構造
4. 画像類似度の計算(Dice係数)
4.1. calc_dice.py
4.2. result_diceindex.csv
5. FSLでやりたい場合
5.1. fsl_dicecalc.sh
1. 目的
2. scikit-learnのインストール
3. データの準備
3.1. label.csv
4. ソースコード
4.1. calc_kappa.py
5. 実行
6. 結果の解釈
7. Neural Network Consoleを使っている場合
7.1. voutput_result.csv
7.2. calc_kappa_nnc.py
1. 目的
2. 準備
2.1. open-visualizationsのダウンロード
2.2. ライブラリのインストール
3. チュートリアル1 (plotnineを用いる場合)
3.1. ライブラリの読み込み
3.2. 保存用のフォルダを用意
3.3. データの読み込み
3.4. データの選択
3.5. プロット
3.6. プロットと直線
3.7. グループごとのプロットの位置を微妙に変える
3.8. プロットの色を変更
3.9. 箱ひげ図 (boxplots)
3.10. バイオリン図 (violin plot)
3.11. 信頼区間 (CI bar)
3.12. 各グループの平均を直線で結ぶ
3.13. プロット・箱ひげ図・バイオリン図・信頼区間
4. チュートリアル2 (matplotlibを使う場合)
4.1. ライブラリの読み込み
4.2. 保存用のフォルダを用意
4.3. データの初期化
4.4. プロット
4.5. プロットと直線
4.6. グループごとのプロットの位置を微妙に変える
4.7. the amount of jitter and create a dataframe containing the jittered x-axis values
4.8. 信頼区間 (CI bar)
4.9. バイオリン図 (violin plot)
4.10. 2群のBeforeとAfterをそれぞれプロット
4.11. さらに信頼区間の追加
4.12. プロット・箱ひげ図・バイオリン図・信頼区間
5. 高画質で保存したい場合
5.1. plotnineの場合
5.2. matplotlibの場合
1. 目的
2. シンプルなヒートマップ
2.1. ライブラリのインポート
2.2. データの読み込み
2.3. 相関行列の計算
2.4. ヒートマップの作成
3. プロットの大きさを相関係数に応じて変える
3.1. heatmapzのインストール
3.2. ライブラリのインポート
3.3. データの読み込み
3.4. ヒートマップの作成
1. 目的
2. 1つの図に3群の結果をプロット
2.1. データ準備
2.2. ソースコード
2.3. 結果確認
2.3.1. Boxplot
2.3.2. Boxplot_with_dot
2.3.3. Violinplot
3. 1つの図に3群の結果を各領域ごとにプロット
3.1. データ準備
3.2. ソースコード
3.3. 結果
4. 1つの図に3つの変数に対して4群の結果を3パターンプロット
4.1. データ準備
4.2. ソースコード
4.3. 結果

1. 目的
2. 準備
2.1. インストール
2.2. 使用データ
2.3. 前処理
3. 神経突起イメージング(NODDI)
3.1. AMICOのセットアップ
3.2. データの読み込み
3.3. 応答関数(response function)の計算
3.4. モデルフィッティング
3.5. 結果
Pythonを使って、AMICOを用いた神経突起イメージング(NODDI)をするために、以下のPythonパッケージをインストールする。
pip3 install dmri-amico
データを次のフォルダ構造で用意する。
Study/
└── Subject
├── DWI.nii.gz # 拡散MRI
├── DWI_mask.nii.gz # 拡散MRIマスク画像
├── bvals # b-values
└── bvecs # b-vectors
NODDI前に、拡散MRIの前処理をする。
Pythonで以下のコマンドを実行。
今回使用するファイル等の変数設定をする。
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()
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 ]
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 ]
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 ]
結果は、「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
画像はこちら。

1. 目的
2. 前準備
3. 予備知識
4. NumPyのインポート
5. データ
6. 重み付けの初期化
7. 学習
8. データの入力
9. 重み(Weight)の勾配計算
10. 損失の計算
11. 重みの更新
12. 実行
12.1. 1_numpy.py
PyTorchのチュートリアルWarm-up: numpyを参考にNumpyを使って、損失(loss)や重み(weight)の計算をする。
PyTorchの特徴の一つである、テンソルとNumpyの違いを理解するための前準備。
PyTorchのインストールはこちらから。
初めて、Google Colaboratoryを使いたい方は、こちらをご覧ください。
脳神経細胞は、樹上突起(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を予測するケースを考えます。
import numpy as np
バッチサイズ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)
乱数を使って重みを初期化します。
# Randomly initialize weights w1 = np.random.randn(D_in, H) w2 = np.random.randn(H, D_out)
学習率を1e-6として、学習回数を500回とします。
learning_rate = 1e-6 for t in range(500):
入力(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)
これより先は、パーセプトロンが予測した値(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)
パーセプトロンが予測した値(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)
計算した勾配(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
以下のコードを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
1. 目的
2. 準備
2.1. DIPYのインストール
2.2. 使用データ
3. 拡散MRIのノイズ除去
3.1. 必要なパッケージをインポート
3.2. 画像およびMPG軸情報の読み込み
3.3. マスク画像の作成
3.4. ギブズのリンギングアーチファクト除去
3.5. NIfTI形式で保存
3.6. 結果
pip3 install dipy
データを次のフォルダ構造で用意する。
Study/
└── Subject
├── DWI.nii.gz # 拡散MRI
├── DWI_mask.nii.gz # 拡散MRIマスク画像
├── bvals # b-values
└── bvecs # b-vectors
Pythonで以下のコマンドを実行。
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
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)
median_otsu関数を用いて、b=0画像からマスク画像を生成する。vol_idxには、b0 volumeのvolume indexを渡す。
maskdata, mask = median_otsu(
data, vol_idx=np.where(bvals == 0)[0])
gibbs_removal関数を用いて、リンギングアーチファクトを除去する。
data_corrected = gibbs_removal(maskdata)
save_nifti関数で、画像をNIfTI形式で保存する。
save_nifti('DWI_degibbs.nii.gz', data_corrected.astype(np.float32), affine)
補正前後の画像は、以下の通り。

最近、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で仮想環境内で構築」です。
概要は以下になります。
venv
2.1 仮想環境の構築
2.2 仮想環境の有効化
2.3 仮想環境の無効化