Numpyと画像とガンマと速度

TL;DR

ファイルに保存された画像データには、ほとんどの場合、ガンマがかかっている。印刷系なら1.8。動画系なら2.2。sRGBならガンマっぽいけれど実はガンマではない、ガンマ2.2に近い処理がされている。

これを元に戻す処理がデガンマ。デガンマしないと、cv2.blurとかで色がおかしくなる。

ガンマのかかったデータがチャンネル当たり8ビットなら、デガンマ後は16ビットで扱う。でないと元のデータが失われ、露骨なトーンジャンプが発生する。

デガンマしたままでデータを引き回すのがリニアワークフロー。データ量は2倍になるが、ガンマ・デガンマは重い処理の上に、ガンマがかかっているかどうか＆ガンマはいくつなのかを間違えることがある。「入り口でデガンマ、出口でガンマ」と決め打ちするメリットは大きい。昔はメモリが足りなかったので、データ量2倍はきつかったのだ。

ガンマの発祥はCRTの特性云々という話があるが、ガンマは人間の知覚に沿った圧縮手法と考えたほうがいい。物理現象（混色）はリニア、知覚はガンマなのだ。

リニアワークフローを貫徹できればいいのだが、時として、16ビットでは困る場合に出くわす。ライブラリが8ビット＆ガンマ前提だったり、速度が欲しかったり。

そのときに愚直に1.8や2.2でガンマする必要があるかどうか考えると、速度がかなり違ってくるかもしれない。

np.powerはnp.sqrtの6倍くらい遅い

%timeit a = np.sqrt(np.arange(1,1000).astype(np.float))
5.89 µs ± 17.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit a = np.power(np.arange(1,1000).astype(np.float), 0.5)
33.1 µs ± 16.9 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

2でガンマすることにして、np.power(x, 0.5)をnp.sqrtに置き換えると、だいぶ速くなる。

追記

人間の知覚に近いガンマは3だが、（2.4弱だった）np.cbrtはなぜかnp.powerより遅い。

%timeit a = np.cbrt(np.arange(1,1000).astype(np.float))
48.9 µs ± 96.3 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

%timeit a = np.power(np.arange(1,1000).astype(np.float), 1/3)
32.1 µs ± 39.8 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

追記2

16ビット→8ビットのLUTは64KBなので、普通のCPUのL1データキャッシュ（Ryzenで32KB）には乗らない。だから平方根を計算するほうが速いだろう、と思い込んでいたが、Pythonではそうでもなかった。

def np_linear16_to_gamma8(img16):
    g = np.zeros(img16.shape, np.uint8)
    less = img16 <= 64
    g[less] = (img16[less] / 8).astype(np.uint8)
    g[~less] = np.sqrt(img16[~less]).astype(np.uint8)
    return g


def cv2_linear16_to_gamma8(img16):
    g = np.zeros(img16.shape, np.uint8)
    less = img16 <= 64
    g[less] = (img16[less] / 8).astype(np.uint8)
    g[~less] = cv2.sqrt(img16[~less].astype(np.float))[:,0].astype(np.uint8)
    return g

_L16_G8_TABLE = np_linear16_to_gamma8(np.arange(65536))

img16 = np.arange(10000,dtype=np.uint16).reshape(100,100)

%timeit a = cv2_linear16_to_gamma8(img16)
69.4 µs ± 259 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

%timeit a = np_linear16_to_gamma8(img16)
93.9 µs ± 356 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

%timeit a=_L16_G8_TABLE[img16]
46.2 µs ± 743 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

np.sqrtが遅い。cv2.sqrtはfloat32とfloat64しか食わないので、uint16からfloat32に型変換するペナルティがあるのに、それでもnp.sqrtより速い。

uint16決め打ちの平方根関数を、コンパイラで最適化が効くように書けば、おそらくLUTよりも速いと思うが、現状のPythonではLUTが一番速いという結論になった。

追記3

デガンマ側も。np.powerよりもcv2.powのほうが6倍くらい速い。

fimg = np.arange(65536).reshape(256, 256) / 65535

%timeit a = cv2.pow(fimg, 1.8)
298 µs ± 1.27 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%timeit a = np.power(fimg, 1.8)
1.84 ms ± 5.83 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

ちなみにcv2.powは1次元のndarrayを食わせると2次元にして返すので要注意。

追記4

なぜcv2.sqrtはこんなに速いのかと思って調べたら、x64のSIMD命令に平方根を求める命令があった。おそらくOpenCVは、Intel IPPを通じてこの命令を使っている。ストリーミング SIMD 拡張命令の算術演算

これではuint16決め打ちの平方根関数など書いても太刀打ちできない。