はじめに
生物系に限らず、画像解析はImageJが普及していますが、最近流行りの機械学習との連携や高度な解析を自動化するにあたって難がありますよね。これからはプログラミングベースでの画像解析を推していきたいです。
Pythonでの画像解析を大幅に効率化したくてこちらのようなライブラリを開発したりしています。せっかくなので、その過程で得た知見などをまとめていこうかなと思います。
今回はサイズの大きいTIFFファイルを効率的に読み込んだり解析したりする方法についてまとめました。
環境
python==3.7.7
dask==2021.6.0
matplotlib==3.3.4
numpy==1.20.3
scikit-image==0.18.1
tifffile==2021.6.14
scikit-imageを使った画像読み込み
画像解析用のライブラリといえばscikit-image。これに含まれているioモジュールのimread関数を使うと、適切なプラグインを選んで画像をnumpy.ndarrayとして読み込んでくれる。tifffileがインストールされていれば簡単にn次元のTIFFファイルを読み込める。
from skimage import io
import matplotlib.pyplot as plt
path = ".../XXX.tif" # パスを格納
img = io.imread(path) # これでOK
img.shape # 形の確認
多くの場合はこれで事足りるが、画像を読み込んで終わりだと少し不自由する場合がある。例えば
- TIFFのメタデータを読み込みたい
- TIFF画像の一部分だけ読み込みたい
などといった場合。今回は2.に着目していく。
TiffFileの使い方
画像を読み込む
モジュールのインポートする。
from tifffile import TiffFile
import numpy as np
よくあるパターンで、withブロック内でファイルを開く。返ってくるファイルハンドルtifはasarrayメソッドで簡単にnumpy.ndarrayに変換できる。
with TiffFile(path) as tif:
img = tif.asarray()
plt.imshow(img) # 画像の表示
画像を保存する
上では何もしていないが、フィルタをかけるなりしたらimwriteを使って画像を保存する。ImageJと互換なメタデータを同時に保存するといった便利なオプションもあるが、今回は割愛する。
from tifffile import imwrite
imwrite(".../save/path/XXX.tif", img) # 保存完了
サイズの大きいTIFFファイルを扱う
メモリマップ
TIFFファイルが大きいと、読み込みに時間がかかってしまったり、そもそもメモリに乗らなかったりして困る。その際に便利なのがメモリマップnumpy.memmapである。詳しくは公式ドキュメントに載っているが、要するにディスク上のファイルをメモリに読み込まない一方でarr[0,3]のようにndarrayを扱っているかのごとく部分アクセスできるもの。
TiffFileは画像をメモリマップとして返してくれるオプションがある。
with TiffFile(path) as tif:
mmap = tif.asarray(out="memmap") # out="memmap"でメモリマップを返してくれる
# 形や型の確認はメモリマップでもできる
print(mmap.shape)
print(mmap.dtype)
メモリマップはいつでもnp.ndarrayに変換できる。
img = np.asarray(mmap) # もちろんすぐに変換してはメモリマップの意味がない！
TIFFを小さくしてから読み込む
さて、大きなファイルが扱いづらいなら、小さくしてから読み込めばいいわけで、メモリマップを使えば簡単に実現できる。例えば、タイムラプス撮影をして、t軸, y軸, x軸の順に(300, 2048, 2048)という形状をしたファイルがあるとする(16 bitなら約2.4 GB)。こういうとき往々にして見たいものが画像のごく一部の場合が多いので、次の手順で読み込む。
- 1枚目だけ読み込んで見たい物の位置を確認する。
- その領域だけ切り出してから読み込む。
まず1枚目を可視化する。
# メモリマップを取得
with TiffFile(path) as tif:
mmap = tif.asarray(out="memmap")
# 1枚目を表示
img0 = np.asarray(mmap[0])
plt.imshow(img0)
仮にy=400～800, x=600～900付近の領域に見たい物があったとしよう。メモリマップでその領域を切り出す(2.4 GB → 72 MB)。
img = np.asarray(mmap[:, 400:800, 600:900])
plt.imshow(img[0]) # 正しく切り出せたかチェック
これでメモリへの読み込みを最小限に抑えられた。
daskの利用して画像を投影する
共焦点顕微鏡のタイムラプス観察など、時間+3次元空間という4次元の画像を扱うとき、z方向に投影して時間+2次元にした方が見やすい場合がある。投影後の画像はサイズが小さくなるので、前述したように何とかしてメモリに全部乗せずに投影後の画像を取得したい。まず思いつく方法は、
-
t=0の画像を読み込む。 - 投影して
listに格納しておく。 -
t=1の画像を読み込む。 - 以下繰り返す。
というもので、確かにメモリに乗せる量を節約できている。しかし、コードが少し煩雑になったりと面倒。加えて、投影だけならまだよいが、他の演算をさせるときにデータサイズが大きいと計算に時間がかかるので、並列処理で効率化したいのが常である。
このときに役立つのがdaskである。
daskとは
daskはnumpy.ndarrayに対応するdask.array.core.Arrayというオブジェクトを持っており、巨大な配列を"chunk"ごとに分けて計算することで、メモリとディスクのデータのやり取りや並列処理を効率的に行ってくる。そして、必要になったときにnumpy.ndarrayやint,floatなどに変換され、値がメモリに残る。daskにはnumpyの関数がほぼ全部実装されており、従来のコードをほぼ変えずに実行できることも特長の一つである。
daskで実際に計算が行われるのは、compute関数が呼び出されたときやnp.array(x)などで明示的にnumpy.ndarrayに変換されるとき、その他if x.mean() > 0: ...のように値が必要になったときである。必要になってから計算を始めることから、英語では「lazyに計算する」(やるべきことを後回しにするというニュアンス？)とかと言われる。
ビッグデータを扱う際のdaskの基本的な流れは
-
daskでできる限りの計算をさせる。 - データが小さくなったところで、
numpy.ndarrayに変換する。
になる。
daskの例
中心極限定理の実験として、「一様乱数100万個の平均値」を1000回計算してヒストグラムを描いてみる。numpyで普通に書くと
data = np.random.random(size=(1000, 1000000)) # 乱数生成
means = np.mean(data, axis=1) # 列ごとに平均値を計算
plt.hist(means, bins=50) # ヒストグラム描画
plt.show()
となるが、dataが8 GBなのでPCによってはアウト。
from dask import array as da
data = da.random.random(size=(1000, 1000000)) # 乱数生成
means = da.mean(data, axis=1) # 列ごとに平均値を計算
plt.hist(means, bins=50) # ヒストグラム描画 (ここでnumpyに変換される)
plt.show()
とすることで、問題なく計算が行われる。npをdaに置き換えただけ。
daskを実際に画像に用いる
da.arrayにメモリマップnp.memmapを渡すと、ノータイムでdask配列に変換される。
with TiffFile(path) as tif:
mmap = tif.asarray(out="memmap")
darr = da.array(mmap)
3次元のタイムラプス撮影では、軸はtzyxの順番に並んでいるので、zで投影するならaxis=1と設定することになる。最大値で投影する場合を例にとると、
img = da.max(darr, axis=0).compute()
でOK。おまけに並列処理でCPUを最大限に活用してくれるので、画像サイズが大きくなるほど効率が良くなる傾向にある。
補足
最近だとdask-imageという、画像解析に特化したdaskのライブラリの開発が進んでいる。今のところ、scipy.ndimageの主要な関数 (gaussian_filter, labelなど) はdask向けに実装されていて、大きな画像の解析が簡単になってきている。
最後に
ほかの拡張子の画像でも、czifileやmrcfileを使えば同様なことができます。
これからも主に画像解析などで色々まとめていきたいと思います。