1. 目的
2. 準備
2.1. data.csv
3. ソースコード
4. 結果
4.1. パラメトリック検定の結果(para_result.csv)
4.2. パラメトリック検定の結果(nonpara_result.csv)
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 仮想環境の無効化
1. 目的
2. 準備
2.1. pydicomのインストール
2.2. dcm2nii / dcm2niixのPATH設定
3. ソースコード
3.1. dcm2niiバージョン
3.2. dcm2niixバージョン
4. 使い方
1. 目的
2. チュートリアル
3. 準備
4. データ取得
5. CSVファイルを読み込み
5.1. Pandasを使ったデータフレーム(DataFrame)のアクセス
5.2. 演算
6. 図示
6.1. Matplotlibを用いたシンプルな方法
6.2. seabornを用いた図示
7. 簡単なデータクリーニングと探索的分析
目的
Pythonとは
Python基本
コメントアウト
一行をコメントアウト
複数をまとめてコメントアウト
文字列の表示:print
四則演算
論理型(Boolean)
論理演算子:and, or
Boolean値の四則演算
比較演算子
文字列処理
Nullオブジェクト:None
変数定義
リスト型
定義
要素を追加
要素を除外
スライス操作
インデックスの確認
結合
ある要素が含まれているか確認
長さ
タプル型
定義
長さ
結合
スライス操作
ある要素が含まれているか確認
辞書型
定義
アクセス
キーの確認
ある要素が含まれているか確認
要素を追加
要素を削除
set型
定義
要素を追加
重複する要素の確認
四則演算
ある要素が含まれているか確認
条件処理:If文
反復処理:For文
反復処理:While文
例外処理
With文
書き込み
読み込み
イテレータ(Iterator)
関数
基本的な使い方
可変長引数:args, kwargs
無名関数:lambda
関数の繰り返し利用:map
条件抽出
内包表記
モジュール(Module)
クラス(Class)
定義
継承
複数の継承
最近、Pythonの勉強をしていますが、いろいろ調べている中で、Matthew Brett氏がMatlabやRに慣れている人のためのPythonの手引き、”Brisk Introduction to Python“を書かれているのを見つけました。既にMatlabやRを知っている方々がPythonのコツを理解するのによいテキストと思います。このため、Brett氏の承諾を得て、研究室のメンバーで翻訳しました。
動的プログラミング言語を知っている人のためのPython入門はこちらからどうぞ。
Ubuntu で Jupyter notebook を入れるとデフォルトでは Firefox が起動します。
Google Chrome にする方法を探しました。
こちらのリンクに方法が記載されていましたが、記載しておきます。
jupyter notebook --generate-config
これで、 ~/.jupyter/jupyter_notebook_config.py が生成されます。
which google-chrome
私の場合は /usr/bin/google-chrome でした。
以下の記載を探します。
# c.NotebookApp.browser = ''
ここを以下に修正します。
c.NotebookApp.browser = '/usr/bin/google-chrome'
これだけです。
2020年2月6・7日に、東京大学駒場キャンパスにおいて、「国際脳MRIプロトコルデータ(HCP data 含む)と精神疾患臨床データの前処理・解析」というテーマにおいて、脳MRI・臨床データ解析チュートリアルが開催されます。ついに日本でHCPデータの解析法などを学べる時が来ました。
主催の東京大学の小池先生から情報をいただきましたので、告知させていただきます。
関心ある方は、こちらをご確認ください。
私は講師としては参加しませんが、HCP readyのLin4Neuroを提供します。
先日、Python 3.x系のAnacondaをLinuxにインストールする必要がありました。
これ自体はなんら問題がなかったのですが、FreeSurferを使っている時にひとつ問題が出ました。
FreeSurferの asegstats2table と aparcstats2table はPython 2.x系で書かれているスクリプトで、Python3系では動きません。
Ubuntuの場合、デフォルトだと
python or python2 –> Python 2.x系
python3 –> Python 3.x系
となっているため、
#!/usr/bin/env python
は、通常はpython2.x系を指します。
しかし、私のように、Python3.x系のAnacondaをいれると
$ python --version
とすると
Python 3.6.5 :: Anaconda, Inc.
となり、python3.x系が呼び出されているのがわかります。
どうしたものかと思いましたが、Anacondaのサイトにきちんと答えが書いてありました。