wide and deep

Light fieldデータセットEPFLをpython3から扱う

python Light field

今回したこと

Light fieldデータセットであるEPFLをpython3から扱い，画像として表示を試みた．
Light field関係はMATLABが優勢っぽいので，pythonで乗り込んでいこうと思う．

本記事で触れたコードはgistに挙げている
gist.github.com

Light field？

Light fieldはとても興味深いものである．
簡単に言うと複数視点から撮影したデータの集合体．
詳しく知りたい人はこの辺を見れば幸せになれると思う．
http://www.cc.kyoto-su.ac.jp/~kano/pdf/paper/2013%20MOC%20Lytro.pdf
netsu-n.mep.titech.ac.jp
btpixel.hatenadiary.org

データセット

Light fieldデータセットは色々存在するが，
主に実世界と仮想空間のデータセットがある．
下記はその一例

他にも下記リポジトリにまとまっているものがある．
github.com

データの取得

今回はLytro Illumで撮影されたEPFLデータセットを扱う．
ここに書いている手順でftpで接続しデータを取得する．
今回取得したデータは下記のもの

EPFL/
┣ 4D_LF/
┃　┣ Landscapes/
┃　┃　┗ Slab_&_Lake.mat
┃　┗ .../
┃
┗ LFR_files/
　　┣ Landscapes/
　　┃　┗ Slab_&_Lake.LFR
　　┗ .../

EPFLデータセットの論文によると

  • LFRファイル...5368x7728のデータ
  • .matファイル...Light Field Matlab Toolbox*1でLFRを加工し，色々なデータを抽出したもの．形状は15x15x434x625x4．15x15のangular resolutionと434x625のspatial resolution，4のカラーチャネル（1つは追加の重み付け成分？）

データをpython3で見ていく

LFRデータ

LFRはimageioというライブラリで読み込めるらしいのでインストール

pip install imageio

取得したSlab_&_Lake.LFRを読み込んでみる

import imageio
import numpy as np

name = 'EPFL/LFR_files/Landscapes/Slab_&_Lake.LFR'
LF = imageio.imread(name)
raw_array = np.array(LF)
h, w = raw_array.shape
print(raw_array.shape)    # (5368, 7728)

読み込んだデータの形状が論文と一致することを確認できた．
このとき各値は0~1.0になっているようだった．
次は表示をしてみる．

import matplotlib.pyplot as plt
plt.figure(figsize=(8, 8))
plt.imshow(raw_array, cmap='gray', vmin = 0, vmax = 1)

f:id:catdance124:20190708221829p:plain
一部を切り出したものはこちら

plt.imshow(raw_array[100:300, 2600:2800], cmap='gray', vmin = 0, vmax = 1)

f:id:catdance124:20190708222156p:plain
小レンズ領域から画像がなっていることがわかる．
しかし，最初から3チャンネルだと思っていたが，そうではないらしい．
調べるとraw画像はデモザイク処理で色付けをするようだ．
下記ページなどを参考に，ベイヤー変換を行う．
colab.research.google.com
www.avaldata.co.jp

データセットからは変換パターンがどれか読み取れなかったので，
下記4パターンを全て試し，GBRGが適当と判断した．
https://www.avaldata.co.jp/solution_imaging/cameralink_tips/images/Bayer_Pattern.png

def simple_demosaic(raw_array, pattern):
    height, width = raw_array.shape
    dms_img = np.zeros((height//2, width//2, 3))
    pattern[pattern == 3] = 1
    dms_img[:, :, pattern[0, 0]] = raw_array[0::2, 0::2]
    dms_img[:, :, pattern[0, 1]] += raw_array[0::2, 1::2]
    dms_img[:, :, pattern[1, 0]] += raw_array[1::2, 0::2]
    dms_img[:, :, pattern[1, 1]] += raw_array[1::2, 1::2]
    dms_img[:, :, 1] /= 2
    return dms_img

#RGGB = np.array([[2, 1], [1, 0]])
#BGGR = np.array([[0, 1], [1, 2]])
#GRBG = np.array([[1, 2], [0, 1]])
GBRG= np.array([[1, 0], [2, 1]])
dms_img = simple_demosaic(raw_array, GBRG)

plt.figure(figsize=(8, 8))
plt.imshow(dms_img,  vmin = 0, vmax = 1)

f:id:catdance124:20190708223449p:plain
4ピクセルから1ピクセルを合成しているので，サイズは元画像の1/4になる．

.matデータ

MATLABから保存されたデータ？の.matはscipyで開けるらしい．
hydrocoast.jp
しかし，バージョン関連のエラーが出た．

from scipy import io
mat= io.loadmat('EPFL/4D_LF/Landscapes/Slab_&_Lake.mat', squeeze_me=True)
# Please use HDF reader for matlab v7.3 files

いつものようにstackoverflowからベストプラクティスを教えてもらい，下記のように実行した．

import h5py
import numpy as np

filepath = 'EPFL/4D_LF/Landscapes/Slab_&_Lake.mat'
mat = {}
f = h5py.File(filepath)
for k, v in f.items():
    mat[k] = np.array(v)

取得したデータの一覧を確認する．

mat.keys()
# dict_keys(['#refs#', 'DecodeOptions', 'GeneratedByInfo', 'LF', 'LFMetadata', 'LensletGridModel', 'RectOptions', 'WhiteImageMetadata'])

この中の'LF'がメインデータ
こっちのLFはなぜか値がよくわからないことになっているので，0~1.0に正規化．
また，扱いやすいように次元軸を入れ替えておく．（角度，角度，空間，空間，色）

LF = mat['LF']
LF.shape    # (4, 625, 434, 15, 15)
LF = LF.transpose(3,4,2,1,0) / LF.max()
LF.shape    # (15, 15, 434, 625, 4)

ここまできたらあとは画像化するのみ
カラーチャネルは4つめの要素が追加の重み付け成分らしいので一旦弾いてRGBとして表示

plt.imshow(LF[7,7,:,:,:3])

f:id:catdance124:20190708225800p:plain

次にLight fieldらしく複数視点の画像を表示してみる．
各視点画像が対応するようにsubplotをする．

plt.figure(figsize=(20,20))
for x in range(15):
    for y in range(15):
        plt.subplot(15, 15, x*15+y+1)
        plt.imshow(LF[x,y,:,:,:3])
plt.show()

f:id:catdance124:20190708230014p:plain
この出力を見ると15x15視点の434x625の画像が取得されていることがよくわかる．
(1,1)や(1,15)の画像が真っ黒なのはケラレってやつだと思う．

.mapデータからremap画像を作成

LFデータには全画素値が入ってるので，そのデータから.LFRデータのような画像を作成してみる．
LFデータは(u,v,x,y,ch)形状なので，(u×x, v×y, 3)の空配列を用意し，格納していくようにする．
出来上がった関数は下記の通り．

def create_remap(LF, inner_n=15):
    x, y = LF.shape[2:4]
    u, v = inner_n, inner_n
    remap = np.zeros((u*x,v*y,3))
    outer_n = (15-inner_n) // 2
    
    for x_i in range(x):
        for y_i in range(y):
            remap[inner_n*x_i:inner_n*(x_i+1), inner_n*y_i:inner_n*(y_i+1)] = \
            LF[outer_n:15-outer_n, outer_n:15-outer_n, x_i, y_i, :3]
    return remap

関数内に存在する，配列コピーの速度を2通りで比較した．
繰り返しを空間ピクセルでおこなうのか，角度ピクセルでおこなうのか，だ．
速度測定はipythonの%%timeでおこなった．

u, v, x, y, _ = LF.shape
remap = np.zeros((u*x,v*y,3))

%%time
for x_i in range(x):
    for y_i in range(y):
        remap[u*x_i:u*(x_i+1), v*y_i:v*(y_i+1)] = LF[:,:,x_i,y_i,:3]
# Wall time: 875 ms

%%time
for u_i in range(u):
    for v_i in range(v):
        remap[u_i::u, v_i::v] = LF[u_i,v_i,:,:,:3]
# Wall time: 3.25 s

結果として，繰り返し回数が多くても（434*625 : 15*15）コピー量が少ない（15*15 : 434*625）ほうが高速となった．
この結果は，繰り返し回数とデータ量により異なると思うが，この関数に置いては上部のコードを採用した．

実際に作成

したものは下記の通り．

remap = create_remap(LF, 15)
print(remap.shape)# (u*x,v*y,3)
plt.imshow(remap)

f:id:catdance124:20190712180132p:plain
拡大したもの

plt.imshow(remap[500:900, 7000:7400])

f:id:catdance124:20190712180244p:plain
また，create_remapでは引数により角度ピクセル(15x15)の内側何ピクセルを使うのかを指定できるようにした．

remap = create_remap(LF, 11)
plt.imshow(remap)

f:id:catdance124:20190712182045p:plain
拡大したもの
マイクロレンズごとのケラレが表示されていないことがわかる

plt.imshow(remap[300:700, 5600:6000])

f:id:catdance124:20190712182114p:plain
create_remap(LF, 1)で作成したもの
普通の画像（中央視点ビュー）が作成できる
f:id:catdance124:20190712183218p:plain
f:id:catdance124:20190712183142p:plain

まとめ

今回はLight fieldデータセットEPFLをpython3から扱う方法についてまとめた．
Light field周辺の情報は日本語だとなかなか落ちていないのでなかなか苦労した．
いずれはHCIデータセットの扱い方についても書きたいと思う．

*1:D. G. Dansereau, “Light-Field Toolbox for Matlab,” December 2015. [Online]. Available: http://www.mathworks.com/matlabcentral/fileexchange/49683-light-field-toolbox-v0-4