サイズが大きいTIFFファイルをPythonで効率的に扱う

# はじめに

生物系に限らず、画像解析はImageJが普及していますが、最近流行りの機械学習との連携や高度な解析を自動化するにあたって難がありますよね。これからはプログラミングベースでの画像解析を推していきたいです。
Pythonでの画像解析を大幅に効率化したくて[こちら](https://github.com/hanjinliu/impy)のようなライブラリを開発したりしています。せっかくなので、その過程で得た知見などをまとめていこうかなと思います。

今回はサイズの大きい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ファイルを読み込める。

```python
from skimage import io
import matplotlib.pyplot as plt

path = ".../XXX.tif" # パスを格納
img = io.imread(path) # これでOK
img.shape # 形の確認
```

多くの場合はこれで事足りるが、画像を読み込んで終わりだと少し不自由する場合がある。例えば

1. TIFFのメタデータを読み込みたい
2. TIFF画像の一部分だけ読み込みたい

などといった場合。今回は2.に着目していく。

# `TiffFile`の使い方

### 画像を読み込む

モジュールのインポートする。

```python
from tifffile import TiffFile
import numpy as np
```

よくあるパターンで、withブロック内でファイルを開く。返ってくるファイルハンドル`tif`は`asarray`メソッドで簡単に`numpy.ndarray`に変換できる。

```python
with TiffFile(path) as tif:
    img = tif.asarray()
plt.imshow(img) # 画像の表示
```

### 画像を保存する

上では何もしていないが、フィルタをかけるなりしたら`imwrite`を使って画像を保存する。ImageJと互換なメタデータを同時に保存するといった便利なオプションもあるが、今回は割愛する。

```python
from tifffile import imwrite
imwrite(".../save/path/XXX.tif", img) # 保存完了
```

# サイズの大きいTIFFファイルを扱う

### メモリマップ

TIFFファイルが大きいと、読み込みに時間がかかってしまったり、そもそもメモリに乗らなかったりして困る。その際に便利なのがメモリマップ`numpy.memmap`である。詳しくは[公式ドキュメント](https://numpy.org/doc/stable/reference/generated/numpy.memmap.html)に載っているが、要するにディスク上のファイルをメモリに読み込まない一方で`arr[0,3]`のように`ndarray`を扱っているかのごとく部分アクセスできるもの。

`TiffFile`は画像をメモリマップとして返してくれるオプションがある。

```python
with TiffFile(path) as tif:
    mmap = tif.asarray(out="memmap") # out="memmap"でメモリマップを返してくれる
# 形や型の確認はメモリマップでもできる
print(mmap.shape)
print(mmap.dtype)
```

メモリマップはいつでも`np.ndarray`に変換できる。

```python
img = np.asarray(mmap) # もちろんすぐに変換してはメモリマップの意味がない！
```

### TIFFを小さくしてから読み込む

さて、大きなファイルが扱いづらいなら、小さくしてから読み込めばいいわけで、メモリマップを使えば簡単に実現できる。例えば、タイムラプス撮影をして、t軸, y軸, x軸の順に`(300, 2048, 2048)`という形状をしたファイルがあるとする(16 bitなら約2.4 GB)。こういうとき往々にして見たいものが画像のごく一部の場合が多いので、次の手順で読み込む。

1. 1枚目だけ読み込んで見たい物の位置を確認する。
2. その領域だけ切り出してから読み込む。

まず1枚目を可視化する。

```python
# メモリマップを取得
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 &rarr; 72 MB)。

```python
img = np.asarray(mmap[:, 400:800, 600:900])
plt.imshow(img[0]) # 正しく切り出せたかチェック
```

これでメモリへの読み込みを最小限に抑えられた。

# `dask`の利用して画像を投影する

共焦点顕微鏡のタイムラプス観察など、時間&plus;3次元空間という4次元の画像を扱うとき、z方向に投影して時間&plus;2次元にした方が見やすい場合がある。投影後の画像はサイズが小さくなるので、前述したように何とかしてメモリに全部乗せずに投影後の画像を取得したい。まず思いつく方法は、

1. `t=0`の画像を読み込む。
2. 投影して`list`に格納しておく。
3. `t=1`の画像を読み込む。
4. 以下繰り返す。
   
というもので、確かにメモリに乗せる量を節約できている。しかし、コードが少し煩雑になったりと面倒。加えて、投影だけならまだよいが、他の演算をさせるときにデータサイズが大きいと計算に時間がかかるので、並列処理で効率化したいのが常である。
このときに役立つのが`dask`である。

### `dask`とは

`dask`は`numpy.ndarray`に対応する`dask.array.core.Array`というオブジェクトを持っており、巨大な配列を"chunk"ごとに分けて計算することで、メモリとディスクのデータのやり取りや並列処理を効率的に行ってくる。そして、必要になったときに`numpy.ndarray`や`int`,`float`などに変換され、値がメモリに残る。`dask`には`numpy`の関数がほぼ全部実装されており、従来のコードをほぼ変えずに実行できることも特長の一つである。

[daskの公式ドキュメント](https://docs.dask.org/en/latest/)

`dask`で実際に計算が行われるのは、`compute`関数が呼び出されたときや`np.array(x)`などで明示的に`numpy.ndarray`に変換されるとき、その他`if x.mean() > 0: ...`のように値が必要になったときである。必要になってから計算を始めることから、英語では「lazyに計算する」(やるべきことを後回しにするというニュアンス？)とかと言われる。

ビッグデータを扱う際の`dask`の基本的な流れは

1. `dask`でできる限りの計算をさせる。
2. データが小さくなったところで、`numpy.ndarray`に変換する。

になる。

### `dask`の例

中心極限定理の実験として、「一様乱数100万個の平均値」を1000回計算してヒストグラムを描いてみる。`numpy`で普通に書くと

```python
data = np.random.random(size=(1000, 1000000)) # 乱数生成
means = np.mean(data, axis=1) # 列ごとに平均値を計算
plt.hist(means, bins=50) # ヒストグラム描画
plt.show()
```

となるが、`data`が8 GBなのでPCによってはアウト。

```python
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`配列に変換される。

```python
with TiffFile(path) as tif:
    mmap = tif.asarray(out="memmap")
darr = da.array(mmap)
```

3次元のタイムラプス撮影では、軸はtzyxの順番に並んでいるので、zで投影するなら`axis=1`と設定することになる。最大値で投影する場合を例にとると、

```python
img = da.max(darr, axis=0).compute()
```

でOK。おまけに並列処理でCPUを最大限に活用してくれるので、画像サイズが大きくなるほど効率が良くなる傾向にある。

###### 補足

最近だと`dask-image`という、画像解析に特化した`dask`のライブラリの開発が進んでいる。今のところ、`scipy.ndimage`の主要な関数 (`gaussian_filter`, `label`など) は`dask`向けに実装されていて、大きな画像の解析が簡単になってきている。

# 最後に

ほかの拡張子の画像でも、`czifile`や`mrcfile`を使えば同様なことができます。
これからも主に画像解析などで色々まとめていきたいと思います。