Python3 & OpenCV で画像処理を学ぶ[5] 〜 AfterEffects/Photoshopにある描画モードを実装する

1. はじめに

　前回は、画像間演算の基礎を学びました。

　今回はその応用として、AfterEffectsなどで利用されている様々な描画モード（合成モード、ブレンドモード）を実装する実験を行い、一部を除き、一通り再現しました。

　注意点として。基本的にインターネット上で挙げられている計算式（例えばこちら https://stackoverflow.com/questions/5919663/how-does-photoshop-blend-two-images-together の回答など）に準じることにしますが、Adobe公式によるものではない（公式は発表していない？）ことも多いので、いくらか誤りを含むかもしれません。また、実装の計算式はアプリケーションごとにも違いがあるため、結果に微妙な違いを含みうることをご了承ください。

　また、描画モード関数の実装は、可読性を最優先した結果いくらか冗長な書き方になっているかと思いますが、ご了承ください。

import numpy as np
from matplotlib import pyplot as plt
import matplotlib.gridspec as gridspec
import seaborn as sns
import cv2

%config InlineBackend.figure_format = 'retina'

2. 合成モード実装実験

2.1 下準備

　あらかじめ、今回描画に使う関数をクラスにまとめておくことにします。

　「二枚の画像を受け取って、ブレンドした画像を返す関数」を渡して生成するクラスです。
合成結果を表示するblend_showメソッドを準備します。

class BlendMode:
    def __init__(self, f):
        """
        合成モードを, ベクトル値関数 f:(x,y)→z の形で与える.
        f は [0~1] かつ同サイズの RGB_img 配列を入出力とすることを要請する.
        """
        self.f = f

    def blend(self, bg_img, fg_img):
        """
        入力 : 前景画像, 背景画像
        出力 : 前景画像 + 背景画像 + 合成結果 の表示
        """
        bg_img = bg_img / 255
        fg_img = fg_img / 255
        
        result = self.f(bg_img, fg_img).clip(0,1)
    
        fig = plt.figure(figsize=(10,9))
        gs = gridspec.GridSpec(3,2)
        
        ax1 = fig.add_subplot(gs[0,0])
        ax2 = fig.add_subplot(gs[0,1])
        ax3 = fig.add_subplot(gs[1:3,:])
        
        sns.set_style('ticks')

        ax1.set_xticks([]), ax1.set_yticks([])
        ax2.set_xticks([]), ax2.set_yticks([])

        ax1.set_title('Background'), ax1.imshow(bg_img)    
        ax2.set_title('Foreground'), ax2.imshow(fg_img)
        ax3.set_title('Result Img'), ax3.imshow(result)
        plt.show()

あとは関数を書いてこのクラスのインスタンスに投げるだけですね。

今回利用する画像素材もまとめて読み込んでおきます。

wiz = cv2.cvtColor(cv2.imread('img_src/Wizdomiot.png'), cv2.COLOR_BGR2RGB)
sozai_grad = cv2.cvtColor(cv2.imread('img_src/RP_grad.png'), cv2.COLOR_BGR2RGB)
sozai_noise = cv2.cvtColor(cv2.imread('img_src/PurpleBlockNoise.png'), cv2.COLOR_BGR2RGB)

imgs_dict = {'Wizdomiot.png':wiz,'sozai_Grad.png':sozai_grad, 'sozai_Noize.png':sozai_noise}

以下に使用素材を示します。

fig = plt.figure(figsize=(18,7))

for (i, img) in enumerate(imgs_dict.items()):
    ax = fig.add_subplot(131 + i)
    ax.set_xticks([]), ax.set_yticks([])
    ax.set_title(str(i) + ". " + img[0])
    ax.imshow(img[1])

plt.show()

また、前回同様、RGB値はすべて 0 と 1 の間の実数として考えます。

$ 0 \leq {RGB} \leq 1$

そして、入出力の画素のRGB値を以下のように書くことにします。

前景（上のレイヤー）:　${RGB}_1$
背景（下のレイヤー）:　${RGB}_2$
出力（演算結果）　　:　${RGB}$

また、${RGB}$ という表記では、実際の計算はR,G,Bチャンネルそれぞれ個別に行われることを示します。例えば、

${RGB} = {RGB}_1 + {RGB}_2$

と書いたとき、実際の計算では、

$R = R_1 + R_2$
$G = G_1 + G_2$
$B = B_1 + B_2$

というように行われるということです。

2.2 「減算」族 〜 重ねると暗くなる

まずは「減算」系の演算です。重ねることで暗くなるという特徴があり、具体的には以下です。

• 比較（暗）
• 乗算
• 焼き込みカラー
• 焼き込みリニア
• カラー比較（暗）

2.2.1 比較(暗) Darken

def Darken(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    # boolean array
    is_BG_darker = bg_img < fg_img
    
    result[is_BG_darker] = bg_img[is_BG_darker]
    result[~is_BG_darker] = fg_img[~is_BG_darker]
    
    return result


darken = BlendMode(Darken)
darken.blend(wiz,sozai_noise)

　比較(暗) Darken は、RGBチャンネルごとに低い方の値を出力にします。

$ {RGB} = min(RGB_1, RGB_2) $

　したがって、上記の例では、背景の魔法陣の白い部分（R,G,B全て強い）と、前景のノイズの明るい部分（Rが強く, G,Bは弱い）とでの比較の結果として、魔法陣のR成分だけが残っている様子などを見てとることができます。

2.2.2 乗算 Multiply

def Multiply(bg_img, fg_img):
    result = bg_img * fg_img
    return result


multiply = BlendMode(Multiply)
multiply.blend(wiz,sozai_noise)

　乗算 multiplyは、二枚の画素値を掛けます。

$ {RGB} = RGB_1 \times RGB_2 $

　RGB値が0~1までであったことを思い出すと、上記の例では、魔法陣の白い部分 ≒ 1 と黒い部分 ≒ 0 が、ノイズの画素値と掛け算されることで、ノイズから魔法陣が浮き上がったかのような結果になっていますね。

2.2.3 焼き込みカラー Color Burn

def ColorBurn(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    bg_inverse = 1 - bg_img
    non_zero = fg_img!=0  # masking array to avoid division by 0
    
    result[non_zero] = 1 - bg_inverse[non_zero]/fg_img[non_zero]
    result[~non_zero] = 0
    
    return result


colorburn = BlendMode(ColorBurn)
colorburn.blend(wiz,sozai_noise)

　焼き込みカラー Color Burnは、「(反転した背景 ÷ 前景) を反転する」という計算を行います。

$ {RGB} = 1- \dfrac{(1 - RGB_2)}{RGB_1} $

　まず、区間[0,1]では、 $(1 - {RGB})$ は反転にあたります。このことは、$y=1-x$ のグラフを描いてみると理解しやすいと思います。
　また、区間[0,1]での割り算を考えるとやや理解しにくいのですが、 シンプルに分数として考えると、以下のようになりそうです。

	分母:小	分母:大
分子:小	≒1	<1
分子:大	>1	≒1

最大値を 1 とすると、分母>分子 のとき以外はすべて1付近になると考えてよさそうです。そこで、「分母のピクセルが明るく、分子のピクセルが暗い」かどうか を基準とすれば多少見通しがよくなりそうです。

　上記の例に即して考えてみましょう。
まず背景を反転すると、「魔法陣が暗くなり、その外側が白っぽくなる」ことは想像がつきます。それを分母=前景のノイズで割ると、「魔法陣の暗さがやや変化し、外側はより明るくなるか100%白に飛ぶ」と考えられます。そこからもう一度反転することで、上記の結果になります。

　結論として、「背景の最も明るい部分は残し、それ以外の部分は急激に暗くなる」描画モードということができそうです。

補足: AfterEffectsの結果と比較したところ、この計算式は、AfterEffects的には「焼き込みカラー（クラシック）」のもののようです。（クラシック）と付く描画モードは基本的に後方互換性のためにあるようですが、「焼き込みカラー（クラシック）」と「焼き込みカラー」では微妙に結果が異なりました。

2.2.4 焼き込みリニア Linear Burn

def LinearBurn(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    bg_inverse = 1 - bg_img
    fg_inverse = 1 - fg_img
    
    result = 1 - (bg_inverse + fg_inverse)
    return result


linearburn = BlendMode(LinearBurn)
linearburn.blend(wiz,sozai_noise)

　焼き込みリニア Linear Burnは、「(反転背景 + 反転前景) を反転する」という計算を行います。

$ {RGB} = 1- \big((1 - {RGB}_1) + (1 - {RGB}_2) \big) $

　反転することで、本来は暗部であったピクセルが明部になります。その状態で加算をして、再び反転して元の世界に戻ってくるということで、「互いの暗部を足し合わせる」というような見方をできそうです。

2.2.5 カラー比較(暗) Darker Color

def DarkerColor(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    if len(bg_img.shape)==3:  # is Input RGB array?
        bg_RGB_sum = bg_img.sum(axis=2, keepdims=True)
        fg_RGB_sum = fg_img.sum(axis=2, keepdims=True)
        
        is_BG_darker = (bg_RGB_sum < fg_RGB_sum).repeat(3,axis=2)  # 配列のサイズ合わせ
        
        result[is_BG_darker] = bg_img[is_BG_darker]
        result[~is_BG_darker] = fg_img[~is_BG_darker]
        
    else:
        result = Darken(bg_img, fg_img)
    
    return result


darkercolor = BlendMode(DarkerColor)
darkercolor.blend(wiz,sozai_noise)

　カラー比較(暗) Darker Color は、RGBの合計値の低い方のピクセルを出力にします。

$ {RGB} = \left\{ \begin{array}{l} {RGB}_1 & (R_1 + G_1 + B_1 < R_2 + G_2 + B_2) \\ {RGB}_2 & (otherwise) \\ \end{array} \right. $

　したがって、「上か下のレイヤーのピクセルがそのまま結果に出る（個々のピクセルは合成で不変）」という特徴があります。

2.3.「加算」族 〜 重ねると明るくなる

• 比較（明）
• スクリーン
• 覆い焼きカラー
• 覆い焼きリニア（加算）
• カラー比較（明）

2.3.1 比較(明) Lighten

def Lighten(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    is_BG_lighter = bg_img > fg_img
    
    result[is_BG_lighter] = bg_img[is_BG_lighter]
    result[~is_BG_lighter] = fg_img[~is_BG_lighter]
    
    return result


lighten = BlendMode(Lighten)
lighten.blend(wiz,sozai_noise)

　比較(明) Lighten は、RGBチャンネルごとに高い方の値を出力にします。

$ {RGB} = max(RGB_1, RGB_2) $

　したがって、上記の例では、背景の魔法陣の白い部分（R,G,B全て強い）がそのまま残っていることと、ノイズにG成分がないために、背景下部の緑色発光のG成分がノイズの上に加わるような結果を見てとることができます。

2.3.2 スクリーン Screen

def Screen(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    result = 1 - ((1 - bg_img) * (1 - fg_img))
    return result


screen = BlendMode(Screen)
screen.blend(wiz,sozai_noise)

　スクリーン Screen は、「反転して乗算し、また反転して戻す」演算です。

$ {RGB} = 1 - (1-RGB_1)*(1-RGB_2) $

　よく加算と比されることが多いですが、計算式の通り、実質的に乗算そのものです。焼き込み(リニア)のように、反転した世界での計算なので、乗算で暗くなる代わりに明るくなっています。

　「加算と比べてスクリーンは白飛びしない」事実は、「スクリーンで白になる = 反転世界の乗算で0になる」と理解するとわかりやすいです。
　加算はすぐ 1 を超えます。しかし掛け算は、少なくともどちらか一方が 0 でない限り 0 にならないため、スクリーン演算ではどちらか一方の画素値が 1 でない限り 1 にはならず、結果として白飛びしません。

2.3.3 覆い焼きカラー Color Dodge

def ColorDodge(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    fg_reverse = 1 - fg_img
    non_zero = fg_reverse!=0
    
    result[non_zero] = bg_img[non_zero]/fg_reverse[non_zero]
    result[~non_zero] = 1
    
    return result

colordodge = BlendMode(ColorDodge)
colordodge.blend(wiz,sozai_noise)

　覆い焼きカラー Color Dodgeは、「背景 ÷ 反転した前景」という計算を行います。

$ {RGB} = \dfrac{(RGB_2)}{1-RGB_1} $

　こちらも除算なので捉えにくいですが、上記の例を眺めながら雰囲気で考えると、「前景のうち明るい部分が反転して低い値になり、それが分母になる」ということと、分数では分母が小さければ値は大きくなることから、「前景の明るい部分が背景に写る」演算とでも理解すればよいでしょうか。
　前景のうち暗い部分は反転して 1 に近い値になり、 1 で割っても値は変わらないので、前景の暗い部分ほど結果に反映されないということのようです。

　また、焼き込みカラーと同様の議論ですが、分数の分母を0に近づけていくと値の増え方が急であるのと同様に、明部から暗部にかけての影響度の変化がきついということが考えられます。実際に、スクリーンなどと比べると、上の例では前景ノイズの暗い部分がほとんど反映されていないことを見て取れるかと思います。

2.3.4 覆い焼きリニアLinear Dodge

def Add(bg_img, fg_img):
    result = bg_img + fg_img
    return result

add = BlendMode(Add)
add.blend(wiz,sozai_noise)

　覆い焼きリニア Linear Dodgeは加算です。1を超えた値は1に抑えられます。

$ {RGB} = {RGB}_2 + {RGB}_1 $

　ただの加算なので、あまり書くことがありません…… このような仰々しい名前が付いているのは、焼き込みリニアと対であることを意識したものかと思われます。

焼き込みリニア :
$ {RGB} = 1- \big((1 - {RGB}_1) + (1 - {RGB}_2) \big) $

2.3.5 カラー比較(明) Lighter Color

def LighterColor(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    if len(bg_img.shape)==3:
        bg_RGB_sum = bg_img.sum(axis=2, keepdims=True)
        fg_RGB_sum = fg_img.sum(axis=2, keepdims=True)
        
        is_BG_lighter = (bg_RGB_sum > fg_RGB_sum).repeat(3,axis=2)
        
        result[is_BG_lighter] = bg_img[is_BG_lighter]
        result[~is_BG_lighter] = fg_img[~is_BG_lighter]
        
    else:
        result = Lighten(bg_img, fg_img)
    
    return result


lightercolor = BlendMode(LighterColor)
lightercolor.blend(wiz,sozai_noise)

　カラー比較(明) Lighter Color は、RGBの合計値の高い方のピクセルを出力にします。

${RGB} = \left\{ \begin{array}{l} {RGB}_1 & (R_1 + G_1 + B_1 > R_2 + G_2 + B_2) \\ {RGB}_2 & (otherwise) \\ \end{array} \right. $

　カラー比較(暗)と全く同様に「上か下のレイヤーのピクセルがそのまま結果に出る（個々のピクセルは合成で不変）」という特徴があります。

2.4「複雑」族 〜 50%グレーを境に、暗部は暗く・明部は明るく

　このカテゴリの描画モードは、基本的に50%グレーを境に暗部をより暗く・明部をより明るくするために、コントラストが上がる傾向にあります。

• オーバーレイ
• ソフトライト
• ハードライト
• リニアライト
• ビビッドライト
• ピンライト
• ハードミックス

　もっと言うと、これまで見てきた描画モードを組み合わせて「明部では加算系の描画モード、暗部では減算系の描画モード」という計算を行うものが多いです。具体的に見ていきましょう。

2.4.1 オーバーレイ Overlay

def Overlay(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    darker = bg_img < 0.5
    bg_inverse = 1 - bg_img
    fg_inverse = 1 - fg_img
    
    result[darker] = bg_img[darker] * fg_img[darker] * 2
    result[~darker] = 1 - bg_inverse[~darker] * fg_inverse[~darker] * 2
    
    return result


overlay = BlendMode(Overlay)
overlay.blend(wiz,sozai_grad)

　オーバーレイ Overlay は、「背景の暗部では乗算、背景の明部ではスクリーン」という演算を行います。

${RGB} = \left\{ \begin{array}{l} 2 ( {RGB}_1 \times {RGB}_2) & ({RGB}_2 < 0.5) \\ 1 - 2 (1 - {RGB}_1 ) ( 1 - {RGB}_2) & (otherwise) \\ \end{array} \right. $

　2をかけているのは、明部と暗部でそれぞれ区間幅が0.5であるための正規化です。

2.4.2 ソフトライト Soft Light

def SoftLight(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    darker = fg_img < 0.5
    
    result[darker] = 2 * fg_img[darker] * bg_img[darker] + 2 * (0.5 - fg_img[darker]) * np.square(bg_img[darker])
    result[~darker] = 2 * (1 - fg_img[~darker]) * bg_img[~darker] + 2 * (fg_img[~darker] - 0.5) * np.sqrt(bg_img[~darker])
    
    return result


softlight = BlendMode(SoftLight)
softlight.blend(wiz,sozai_grad)

　ソフトライト Soft Light は若干ややこしく、またアプリケーションごとに微妙に異なるようなのですが、Photoshopの場合は、「前景の画素値に応じて、背景にチャンネル単位のガンマ補正を行う」という処理が行われているようです。

${RGB} = \left\{ \begin{array}{l} 2 \times {RGB}_1 \times {RGB}_2 + 2 \times (0.5-{RGB}_1) \times {RGB}_2^2 & ({RGB}_1 < 0.5) \\ 2 \times (1 - {RGB}_1) \times {RGB}_2 + 2 \times ({RGB}_1-0.5) \times {RGB}_2^\frac{1}{2} & (otherwise) \\ \end{array} \right. $

まず、$ RGB_2^{2} $, $ RGB_2^\frac{1}{2} $ は、それぞれガンマ補正とみなすことができます。

ガンマ補正とは、(出力) = ${RGB}^\gamma$ とする変換で、$\gamma < 1$ のときは階調が全体的に明るくなり、$\gamma > 1$ のときは全体的に暗くなるものでした。（$0\leq{RGB}\leq1$ であるため）

それを踏まえて、上の式を日本語で直感的に書くと、

背景 $ RGB_2 $ に対するガンマ補正 ${RGB}^\gamma$を考える。
${RGB}_1 \leq 0.5$ のとき、$\gamma = 2$ (暗) と $\gamma = 1$ の間を、${RGB}_1$の値で内分する。
${RGB}_1 \geq 0.5$ のとき、$\gamma = \frac{1}{2}$ (明) と$\gamma = 1$ の間を、${RGB}_1$の値で内分する。
ただし、${RGB}_1 = 0.5$ のときには ${RGB} = {RGB}_2$ で連続的に繋がるように、内分の比をそれぞれ $0.5$ ずらして調整する。

ということになると思います。

そう考えると、この描画モードの特徴は、

最大でも$RGB_2^\frac{1}{2}$まで、最低でも$RGB_2^2$までに値が収まる。

ということが言えそうです。したがって、他の描画モードと比べると穏やかな結果になるということですね。

2.4.3 ハードライト Hard Light

def HardLight(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    darker = fg_img < 0.5
    bg_inverse = 1 - bg_img
    fg_inverse = 1 - fg_img
    
    result[darker] = 2 * bg_img[darker] * fg_img[darker]
    result[~darker] = 1 - 2 * bg_inverse[~darker] * fg_inverse[~darker]
    
    return result


hadlight = BlendMode(HardLight)
hadlight.blend(wiz,sozai_grad)

　ハードライト Hard Light は、オーバーレイで前景と背景を入れ替えた計算を行います。つまり、オーバーレイは背景の明暗で場合分けをしていましたが、ハードライトでは前景の明暗で場合分けで行うということです。

${RGB} = \left\{ \begin{array}{l} 2 ( {RGB}_1 \times {RGB}_2) & ({RGB}_1 < 0.5) \\ 1 - 2 (1 - {RGB}_1 ) ( 1 - {RGB}_2) & (otherwise) \\ \end{array} \right. $

というわけで、上の例でレイヤーを入れ替えてみると、

hadlight.blend(sozai_grad, wiz)

先ほどのオーバーレイと同じ結果が出てきました。

2.4.4 ビビッドライト Vivid Light

def VividLight(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    darker = fg_img < 0.5
    
    result[darker] = ColorBurn(bg_img[darker], 2*fg_img[darker])
    result[~darker] = ColorDodge(bg_img[~darker], 2*(fg_img[~darker]-0.5))
    
    return result


vividlight = BlendMode(VividLight)
vividlight.blend(wiz,sozai_noise)

　ビビッドライト VividLight は、「前景の暗部では焼き込みカラー、前景の明部では覆い焼きカラー」という演算を行います。

${RGB} = \left\{ \begin{array}{l} 1-\dfrac{(1-RGB_2)}{2 \times {RGB}_1} & ({RGB}_1 < 0.5) \\ \dfrac{RGB_2}{\big(1-2 ( RGB_1 -0.5 ) \big)} & (otherwise) \\ \end{array} \right. $

「除算の性質上、暗部だけ強く乗る焼き込みカラー/明部だけ強く残る覆い焼きカラー」の組み合わせによって、極めてコントラストが高くなります。

2.4.5 リニアライト Linear Light

def LinearLight(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    darker = fg_img < 0.5
    
    result[darker] = LinearBurn(bg_img[darker], 2*fg_img[darker])
    result[~darker] = Add(bg_img[~darker], 2*(fg_img[~darker]-0.5))
    
    return result


linearlight = BlendMode(LinearLight)
linearlight.blend(wiz,sozai_noise)

　リニアライト Linear Light は、「前景の暗部では焼き込みリニア、前景の明部では覆い焼きリニア」という演算を行います。

${RGB} = \left\{ \begin{array}{l} 1-\big( (1-RGB_2) + (1 - 2 \times {RGB}_1) \big) & ({RGB}_1 < 0.5) \\ RGB_2 + 2( RGB_1 -0.5 ) & (otherwise) \\ \end{array} \right. $

焼き込みリニアは「反転して加算」、覆い焼きリニアはそのまま加算でした。 前景の明部と暗部が、背景にそれぞれ足しこまれるという感じですね。

2.4.6 ピンライト Pin Light

def PinLight(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    darker = fg_img < 0.5
    
    result[darker] = Darken(bg_img[darker], 2*fg_img[darker])
    result[~darker] = Lighten(bg_img[~darker], 2*(fg_img[~darker]-0.5))
    
    return result


pinlight = BlendMode(PinLight)
pinlight.blend(wiz,sozai_noise)

　ピンライト Pin Light は、「前景の暗部では比較(暗)、前景の明部では比較(明)」という演算を行います。

${RGB} = \left\{ \begin{array}{l} min( RGB_2, 2 \times {RGB}_1 ) & ({RGB}_1 < 0.5) \\ max\big( RGB_2, 2( RGB_1 -0.5 )\big) & (otherwise) \\ \end{array} \right. $

前景または背景の値が残るため、どぎつくはなりにくい感じでしょうか。。。

2.4.7 ハードミックス Hard Mix

def HardMix(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    vivid = VividLight(bg_img, fg_img)
    
    darker = vivid < 0.5
    
    result[darker] = 0
    result[~darker] = 1
    
    return result


hardmix = BlendMode(HardMix)
hardmix.blend(wiz,sozai_noise)

　ハードミックス HardMix は、ビビッドライトの演算結果をチャンネルごとに二値化します。

${RGB} = \left\{ \begin{array}{l} 0 & \big( VividLight({RGB}_2, {RGB}_1) < 0.5\big) \\ 1 & (otherwise) \\ \end{array} \right. $

結果、最大でも8色になることが著しい特徴です。

2.5. 「その他」族 〜 差を使ったり、その他

その他です。あまり普段は使わないものが多いです。

• 差
• 除外
• 減算
• 除算

2.5.1 差 Difference

def Difference(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    result = np.absolute(bg_img - fg_img)
    
    return result


difference = BlendMode(Difference)
difference.blend(wiz,sozai_noise)

差 Difference は、前景と背景の値の差をとります。

$ {RGB} = abs({RGB}_2 - {RGB}_1) $

absは絶対値のことです。差ですね……距離と言い換えてもいいかもしれません。（マンハッタン距離ですね）

2枚の画像が完全に同じであるかを、差で重ねて完全に黒くなるかどうかで判定するという使い方があります。

2.5.2 除外 Exclusion

def Exclusion(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    result = bg_img + fg_img - 2 * bg_img * fg_img
    
    return result


exclusion = BlendMode(Exclusion)
exclusion.blend(wiz,sozai_noise)

除外 Difference は、前景と背景の相加平均と相乗平均の差をとって2倍するような演算をします。

$ {RGB} = 2 \times \big( \dfrac{{RGB}_2 + {RGB}_1}{2} - {RGB}_2 \times {RGB}_1 \big) $

参考までに、相加平均・相乗平均について以下に記します。

相加平均: $\dfrac{A+B}{2}$

相乗平均: $\sqrt{AB}$

また、

大小関係: $\dfrac{A+B}{2} \geq \sqrt{AB}$ （等しくなるのは $A=B$ のとき）

であり、「相加平均 - 相乗平均」は$A=B$のときに$0$となるので、相加・相乗の差は、 $A$ と $B$ との"距離"を別の形で表現したものと見ることができそうです。
（「相加 - 相乗」は三角不等式を満たさないので、厳密には 「距離」 と呼べませんが……） 上の式では√がついてないので、実際には相乗"平均"でもないのですが……それはそうとして、上の結果を先ほどの「差」と比較して見てみると、こちらも二枚のRGB値の間の "距離" が、また別の形で結果に現れているとみることはできそうです。

相乗平均でないために $A=B$ のときに $0$ となってくれるわけではないので、画像の同一性判定には使えません。

2.5.3 減算 Substract

def Substract(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    result = bg_img - fg_img
    result[result<0] =  0
    
    return result


substract = BlendMode(Substract)
substract.blend(wiz,sozai_noise)

減算 Substract は、背景から前景を引き算します。

$ {RGB} = {RGB}_2 - {RGB}_1 $
（ただし、0を下回った値は0）

ノイズが赤っぽいので、赤成分が引かれて緑っぽくなっています。

2.5.4 除算 Division

def Division(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    non_zero = fg_img!=0
    
    result[non_zero] = bg_img[non_zero]/fg_img[non_zero]
    result[~non_zero] = 1
    
    return result


division = BlendMode(Division)
division.blend(wiz,sozai_noise)

除算 Division は、前景で背景を割ります。

$ {RGB} = \dfrac{{RGB}_2}{{RGB}_1} $

割り算は分子・分母によって値が急激に変わったり、あっさり1を超えたりするので扱いが難しいですね……。

2.6 「HSL」族 （失敗）

このカテゴリでは、RGBではなく 色相/彩度/輝度 に関わる合成を行うのですが、輝度の扱いがHSVのVなどとは異なる（人間の分光比視感度を考慮した）もののようで、色空間の変換がすぐにはできず、今回再現に成功しませんでした。

• 色相
• 彩度
• カラー
• 輝度

一応参考までに、書いたコードと結果および、AfterEffectsでの結果を併記しながら、以下に見ていきます。

2.6.1 色相 Hue

def Hue(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    fg_img = (fg_img*255).astype('uint8')
    bg_img = (bg_img*255).astype('uint8')
    
    bg_hsv = cv2.cvtColor(bg_img, cv2.COLOR_RGB2HSV)
    fg_hsv = cv2.cvtColor(fg_img, cv2.COLOR_RGB2HSV)
    
    bg_hsv[:,:,0] = fg_hsv[:,:,0]
    
    img_rgb = cv2.cvtColor(bg_hsv, cv2.COLOR_HSV2RGB)
    
    result = img_rgb/255
    
    return result

hue = BlendMode(Hue)
hue.blend(wiz,sozai_noise)

下のレイヤーの輝度と彩度を維持したまま、上のレイヤーの色相だけそのまま移します。

HSVでは色相ごとの輝度の違いは考慮されていないので、HSV系のエフェクトで色相を動かすよりこちらの方が良さそうな感じがします。

AfterEffectsでの結果 :

2.6.2 彩度 Saturation

def Saturation(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    fg_img = (fg_img*255).astype('uint8')
    bg_img = (bg_img*255).astype('uint8')

    bg_hsv = cv2.cvtColor(bg_img, cv2.COLOR_RGB2HSV)
    fg_hsv = cv2.cvtColor(fg_img, cv2.COLOR_RGB2HSV)
    
    bg_hsv[:,:,1] = fg_hsv[:,:,1]
    
    img_rgb = cv2.cvtColor(bg_hsv, cv2.COLOR_HSV2RGB)
    
    result = img_rgb/255
    return result

saturation = BlendMode(Saturation)
saturation.blend(wiz,sozai_noise)

下のレイヤーの輝度と色相を維持したまま、上のレイヤーの彩度だけそのまま移します。

実際的には、情報量を増やすために彩度にムラをつけるというような用途になるでしょうか。

上の結果を見ると、100%白の場合は本来彩度がないために？ややバグっているように見えます。

AfterEffectsでの結果 :

2.6.3 カラー Color

def Saturation(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    fg_img = (fg_img*255).astype('uint8')
    bg_img = (bg_img*255).astype('uint8')
    
    
    bg_hsv = cv2.cvtColor(bg_img, cv2.COLOR_RGB2HSV_FULL)
    fg_hsv = cv2.cvtColor(fg_img, cv2.COLOR_RGB2HSV_FULL)
    
    bg_hsv[:,:,:2] = fg_hsv[:,:,:2]
    
    img_rgb = cv2.cvtColor(bg_hsv, cv2.COLOR_HSV2RGB_FULL)
    
    result = img_rgb/255
    return result

saturation = BlendMode(Saturation)
saturation.blend(wiz,sozai_noise)

下のレイヤーの輝度を維持したまま、上のレイヤーの色相と彩度をそのまま移します。

イラストなどを描くときに使えそうな感じでしょうか。

AfterEffectsでの結果 :

2.6.4 輝度 luminosity

def Saturation(bg_img, fg_img):
    result = np.zeros(bg_img.shape)
    
    fg_img = (fg_img*255).astype('uint8')
    bg_img = (bg_img*255).astype('uint8')
    
    
    bg_hsv = cv2.cvtColor(bg_img, cv2.COLOR_RGB2HSV_FULL)
    fg_hsv = cv2.cvtColor(fg_img, cv2.COLOR_RGB2HSV_FULL)
    
    bg_hsv[:,:,2] = fg_hsv[:,:,2]
    
    img_rgb = cv2.cvtColor(bg_hsv, cv2.COLOR_HSV2RGB_FULL)
    
    result = img_rgb/255
    return result

saturation = BlendMode(Saturation)
saturation.blend(wiz, sozai_noise)

下のレイヤーの色相・彩度を維持したまま、上のレイヤーの輝度をそのまま移します。

あまり使いどころがなさそうな感じがしますが、上下を入れ替えると「カラー」と同じになります。

AfterEffectsでの結果 :

3. おわりに

3.1 まとめ

今回は、自分で実装することを通じて、普段利用している描画モードについて理解を深めることに成功しました。 また、ブログという形でこうして書くことで、各描画モードをそれぞれ個別に考察する機会ができたのがよかったです。

また、描画モードは数が多く見えますが、処理的には対応関係にあるものが多いということがわかりました。 そこで、描画モード間の関係を以下のような表に整理してみましたので、参考までに。

3.2 課題

色相/彩度/輝度について調べ直したところ、RGBからの変換方式によってHSLやHSVなどと色々異なるモデルが存在することがわかりました。 今回はもう体力がないので、このあたりの考察は今後の課題とします。

https://en.wikipedia.org/wiki/HSL_and_HSV#Hue_and_chroma

4. 参考文献

• https://stackoverflow.com/questions/5919663/how-does-photoshop-blend-two-images-together
• http://fe0km.blog.fc2.com/blog-entry-77.html
• https://dunnbypaul.net/blends/
• https://en.wikipedia.org/wiki/Blend_modes
• https://ofo.jp/osakana/cgtips/blendmode.phtml
• https://handywebdesign.net/2012/07/photoshop-blend-mode-part1/

Python3 & OpenCV で画像処理を学ぶ[4] 〜 アルファブレンドとエンボス効果（画像間演算の基礎）

1. はじめに

　前回は、一つの画像の画素値に対して個々に変換を行うトーンカーブについて学びました。

　今回は、二つの画像の画素値を受け取って変換を行い、一つの画素値として返す処理を学びます。具体的には、画像の上に画像を重ねた時の処理についてです。
　その代表的な例として、アルファブレンディングやエンボス効果などについてまとめ、実装を行ってみたいと思います。

import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns
import cv2

%config InlineBackend.figure_format = 'retina'

2. 画像ブレンドの基礎

2.1. 画像間演算

　二枚(以上)の画像が与えられたとき、同じ位置にある各画像の画素値に対して何らかの演算を行い、一つの画素値として出力することを画像間演算と呼びます。そのような処理にはどのようなものがあるでしょうか。

　画面上重なる位置にある（＝入力される）2つの画素の値を $x, y$と、 出力される画素値を $z$ とおきます。シンプルな例では、

• 加算
  $z = x + y$

• 乗算
  $z = xy$

• 単純平均
  $z = \dfrac{x + y}{2}$

などがありえそうです。

「$x,y$ の二つの入力が与えられることで、 $z$ の値が一通りに決まる」関係ということになるので、二変数関数 $z=f(x,y)$ として考えられそうです。

　以下、より具体的に見ていくことにします。

2.2. アルファブレンディング

　(R,G,B) の3つのチャンネルの他、不透明度を示すアルファ値 A を加えて (R, G, B, A)とすることがあり、Photoshopなどのレイヤー構造を持つソフトウェアでは、アルファ値によって下のレイヤーが透けたりすることは周知だと思います。

　アルファ値による透過に見えるこの処理（アルファブレンディング）が、画像間演算としてはどのようなものになるのかを見ていくことにします。

　以下では簡単のため、RGB値とアルファ値はすべて 0 と 1 の間の実数として考えます。単に255で割れば[0,255]と対応します。

$ 0 \leq {RGB} \leq 1$
$0 \leq \alpha \leq 1 $

そして、入出力の画素のRGB値, アルファの値を以下のように書くことにします。

前景（上のレイヤー）:　${RGB}_1$, $\alpha_1$
背景（下のレイヤー）:　${RGB}_2$, $\alpha_2$
出力（演算結果）　　:　${RGB}$, $\alpha$

2.2.1 $\alpha_2=1$ と仮定する場合

　下のレイヤーの画素が完全に不透明（$\alpha_2=1$）なときの式は、以下のようになります。

下のレイヤーが完全に不透明なため、上にどんな不透明度の画素が乗っても当然出力は $\alpha=1$ です。

　また、RGB値については、幾何的に${RGB}_1,{RGB}_2$をRGB空間中の二点と考えたとき、それらを $\alpha_1 : (1-\alpha_1)$ で内分した点が $ {RGB} $ になると考える事ができます。

以下の図のようなイメージです。

つまり、上のレイヤーの不透明度 $\alpha_1$ が大きければ、出力が上のレイヤーの ${RGB}_1$ に近くなり、逆に$\alpha_1$ が小さければ、${RGB}_2$ に近くなるということですね。

この演算によって、「上のレイヤーが透けている」ような見え方が実現できている、ということです。

2.2.2 $\alpha_2=1$ とは限らない一般の場合

次に、下のレイヤーが半透明（$\alpha_2\neq1$）なときの式は、以下のようになります。

（ただし, $\alpha=0$ つまり出力画素が完全な透明のときは ${RGB}=0$ とする.）

ここで、

$\alpha = \ \alpha_1 + \alpha_2 (1 - \alpha_1)$

の式を理解するために、この式を

$z = x + y (1 - x)$

という二変数関数であるとみなして、三次元作図をしてみることにします。

# jupyterで三次元プロットを動的に動かすには、以下のマジックコマンドでグラフを外部出力する必要あり
%matplotlib 

N= 20
x = np.linspace(0, 1, N)
y = np.linspace(0, 1, N)
X,Y = np.meshgrid(x,y)
P = np.c_[np.ravel(X), np.ravel(Y)]

def alpha(A):
    Z = A[0] + A[1] * (1-A[0])
    return Z

Z = np.apply_along_axis(alpha, 1, P).reshape(X.shape)

sns.set_style('ticks')

fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
surf = ax.plot_surface(X,Y,Z,cmap='binary',linewidth=4)
ax.set_xlabel('foreground alpha')
ax.set_ylabel('background alpha')
ax.set_zlabel('blended')
plt.show()

すると、以上のようになりました。
左奥方向が背景のアルファ値$\alpha_2$、右方向が前景のアルファ値$\alpha_1$、上方向がブレンドした結果の$\alpha$ です。

　こうしてみると、 $\alpha = \ \alpha_1 + \alpha_2 (1 - \alpha_1)$ という式は、「半透明な2枚の層の重ね合わせ」に対して自然に要請される性質、すなわち「常に$max(\alpha_1,\alpha_2)\leq \alpha$ を満たし(合成で減らない)」、かつ「少なくとも$\alpha_1,\alpha_2$のどちらか一方が $1$ でない限り結果も$\alpha=1$にならない」という条件を満たすような演算であることが視覚的にわかります。

2.2.3 プリマルチプライド

ところで、先ほどの式

$ \left\{ \begin{array}{l} \alpha &= \alpha_1 + \alpha_2 (1 - \alpha_1) \\ {RGB} &= \big({RGB}_1 \alpha_1 + {RGB}_2 \alpha_2 (1 - \alpha_1) \big) \div \alpha \\ \end{array} \right. $

のRGBの行では、${RGB}_1 \alpha_1 + {RGB}_2 \alpha_2$ というようにRGB値にその画素のアルファ値を乗算しています(→透明な部分が暗くなる)。

　これをアルファブレンドの際に行わなくて良いように、あらかじめ画像のRGB値をアルファ値を乗算したものにしておくという方式があります（乗算済みアルファ、プリマルチプライド）。

そのようにした場合、アルファブレンドの計算式は以下のようになります。

$ \left\{ \begin{array}{l} \alpha &= \ \alpha_1 + \alpha_2 (1 - \alpha_1) \\ {RGB} &= \ {RGB}_1 + {RGB}_2 (1 - \alpha_1) \\ \end{array} \right. $

かなり式が簡潔になりました。合成のたびにアルファ値を乗算したり割ったりする手間が省けるため、計算の高速化が見込めそうです。

2.2.4 アルファブレンド実装実験

　理論的な話が多くなったので、ここで手を動かす事にします。上記の式を元に、アルファブレンドを実装をしてみます。

以下は、それぞれ下半分・右半分にマスクを書けた二枚の画像をプリマルチプライしたのち、ブレンドを行うコードです。

%matplotlib inline
"""
二枚の画像のアルファ値を調整し下半分・右半分をマスキング。
プリマルチプライしたのちにアルファブレンドを実行して表示。
"""

# imread(path,-1)でBGRAとして読み込み. ともに(720, 1280, 4)の配列
pdr = cv2.cvtColor(cv2.imread('img_src/pandra.png',-1),cv2.COLOR_BGRA2RGBA)
sb = cv2.cvtColor(cv2.imread('img_src/SCARY BANQUET1.png',-1),cv2.COLOR_BGRA2RGBA)

# [0,1] に正規化
pdr = pdr / 255
sb = sb / 255

y = pdr.shape[0]
x = pdr.shape[1]

# やや強引だが、pdrを画面の上1/3から2/3にかけてα=1→0となるようにする。sbは左から右で
# pdrを背景, sbを前景にする想定
pdr[:,:,3] = np.linspace(1.5,-0.5,y).repeat(x).reshape(y,x).clip(0,1)
sb[:,:,3] = np.linspace(1.5,-0.5,x).repeat(y).reshape(x,y).T.clip(0,1)

# Premultiply(今回は一回だけの合成なのであまり意味がないが)
for (i,ch) in enumerate(['R','G','B']):
    pdr[:,:,i] = pdr[:,:,i] * pdr[:,:,3]
    sb[:,:,i] = sb[:,:,i] * sb[:,:,3]

blended = np.zeros(pdr.shape)

# alphaの式
blended[:,:,3] = sb[:,:,3] + pdr[:,:,3] * (1 - sb[:,:,3])

# RGBの式
for (i,ch) in enumerate(['R','G','B']):
    blended[:,:,i] = sb[:,:,i] + pdr[:,:,i] * (1 - sb[:,:,3])



gs = gridspec.GridSpec(3,2)
fig = plt.figure(figsize=(10,9))
ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[0,1])
ax3 = fig.add_subplot(gs[1:3,:])

ax1.set_xticks([]), ax1.set_yticks([])
ax2.set_xticks([]), ax2.set_yticks([])

# アルファは無視して表示
ax1.set_title('Background')
ax1.imshow(pdr[:,:,:3])

ax2.set_title('Foreground')
ax2.imshow(sb[:,:,:3])

ax3.set_title('Blended')
ax3.imshow(blended[:,:,:3])

plt.show()

以上のようになりました。コードの通り、上に表示された画像自体はアルファチャンネルを持っていません。
つまり、不透明度はアルファ値によって直接表現されているわけではなく、アルファ値をRGB値に乗算することによって表現されるということです。

次に、上の画像のプリマルチプライドを解除してみます。つまり、RGB値をアルファ値で割ります。

# プリマルチプライドの解除つまり除算をするので、ゼロ除算回避用のマスク
mask = blended[:,:,3]!=0

upmply = np.zeros(blended.shape)

for (i,ch) in enumerate(['R','G','B']):
    upmply[mask,i] = blended[mask,i] / blended[mask,3]

upmply = (upmply*255).astype('uint8')

plt.figure(figsize=(16,9))
plt.title('un-premultiplied')
plt.imshow(upmply[:,:,:3])
plt.show()

以上のようになりました。上の画像（本来のRGB値）にアルファ値が乗算されることで、半透明に見えていたというわけですね。
また、二枚の画像が重なる部分では、先ほど説明した「${RGB}_1,{RGB}_2$ を $\alpha_1 : (1-\alpha_1)$ で内分した点が $ {RGB} $ になる」という事の意味を視覚的に確かめることができます。

次に、アルファ値の計算式の確認として、 α=1 となるピクセルのみを表示することにしてみます。

partialy_opac = blended[:,:,3]!=1

opac = blended

for (i,ch) in enumerate(['R','G','B']):
    opac[partialy_opac,i] = blended[partialy_opac,i] * 0

opac = (opac*255).astype('uint8')
    
plt.figure(figsize=(16,9))
plt.title('only α=1 pixels')
plt.imshow(opac[:,:,:3])
plt.show()

以上のようになりました。二枚の画像のどちらか一方がα=1となる領域だけが残ることがわかります。

この結果から、「少なくとも$\alpha_1,\alpha_2$のどちらか一方が $1$ でない限り結果も$\alpha=1$にならない」という事実を確かめることができました。

2.2.5 アルファブレンドのまとめ

　結論として、アルファ値というのは見かけ上「不透明度」であり、そのように捉えてほぼ問題ないのですが、計算上はあくまで「画像間演算の際、それぞれのRGB値に対して乗算される係数」に過ぎず、不透明度をそれ自身で直接表現しているものではないということです。

　以上のことを原理的に理解すれば、例えば「不透明度50%のレイヤーを二枚重ねても不透明度100%にならない」などの個別具体的な現象の意味も了解できるかと思います。

2.3. エンボス効果

　次はエンボス効果です。こちらはアルファブレンドに比べると話題としてはややおまけ程度なのですが、画像間演算を別の側面から眺める例として。

　画像を浮き上がったように見せるエンボス効果は、以下の手続きで実現する事ができます。

「ある元画像 $A$ と、$A$の明度を反転して少し斜めに移動した$A'$を加算合成し、$0.5$を引く。ただし、$0$を下回ったり$1$を超えた値は、$[0,1]$の範囲にクリップする」

　早速、具体的に実装していきたいと思います。 まずは元画像です。

# グレースケールにしておく
img = cv2.cvtColor(cv2.imread('img_src/lyrith.png',0),cv2.COLOR_GRAY2RGB)
img = img / 255

plt.figure(figsize=(16,9))
plt.imshow(img)
plt.show()

　次に、この画像の明度を反転します。前回の記事のように、シンプルに階調変換関数を作って適用することにします。

def invert(pixel):
    return 1 - pixel

img_2 = np.array([invert(a) for a in img])

plt.figure(figsize=(16,9))
plt.imshow(img_2)
plt.show()

　反転できました。

　次は斜めへの移動です。

　あるピクセルの座標 $\mathbf{P} = \left( \begin{array}{c} p_x \\ p_y \\ \end{array} \right)$ を、 $ \left( \begin{array}{c} x \\ y \\ \end{array} \right)$方向へ $ \left( \begin{array}{c} \Delta x \\ \Delta y \\ \end{array} \right)$ だけ移動させるには、以下の変換行列 $\mathbf{T}$ を用います。

$ \mathbf{T} = \left( \begin{array}{ccc} 1 & 0 & \Delta x \\ 0 & 1 & \Delta y \\ \end{array} \right) $

　すなわち、$\mathbf{P}$ に 変換行列 $\mathbf{T}$ を左からかけることで、移動後の座標$\mathbf{P'}$ を

$ \begin{align} \mathbf{P'} &= \mathbf{TP} \\ &= \left( \begin{array}{ccc} 1 & 0 & \Delta x \\ 0 & 1 & \Delta y \\ \end{array} \right) \left( \begin{array}{c} p_x \\ p_y \\ \end{array} \right) \\ &= \left( \begin{array}{c} p_x + \Delta x \\ p_y + \Delta y \\ \end{array} \right) \end{align} $

として、得ます。

　OpenCVでこれを行うには、cv2.warpAffine(img, Matrix,(colnum,rownum)) 関数を用います。 第二引数のMatrixに、変換行列をfloat32のnumpy配列として与えることで、変換後の画像を得ることができます。

# 移動量
dx = -5
dy = -5

# 変換行列
M = np.float32([[1,0,dx],[0,1,dy]])

rownum,colnum = img_2.shape[:2]

plt.figure(figsize=(16,9))
img_2 = cv2.warpAffine(img_2, M,(colnum,rownum)) # (列数, 行数)の順であることに注意
plt.imshow(img_2)
plt.show()

　平行移動ができました。(移動が微小なのでわかりにくいですが、右下を見ると少しずれています)

　最後に、これと元画像を重ねて $0.5$ を引きます。

add = img + img_2
emboss = add - 1/2

emboss = emboss.clip(0,1)

plt.figure(figsize=(16,9))
plt.imshow(emboss)
plt.show()

　エンボス効果ができました。リリスだと映えますね。

　このような処理でエンボス効果になる仕組みは、以下の図のようなイメージです。

元画像とそれを反転した画像を単に加算合成すると何も見えなくなりますが、少しの平行移動によって凹凸だけが残るようになる、という感じですね。

3. おわりに

　今回は、アルファブレンドとエンボス効果について学びながら実装することを通して、画像ブレンド（画像間演算）の基礎を理解することができました。

　次回↓は、（やや寄り道気味ですが、）Photoshopなどにある描画モードを一通り実装します。

optie.hatenablog.com

4. 参考文献

• 「ディジタル画像処理」の5章
• アルファブレンド https://ja.wikipedia.org/wiki/%E3%82%A2%E3%83%AB%E3%83%95%E3%82%A1%E3%83%96%E3%83%AC%E3%83%B3%E3%83%89
• アルファチャンネル https://ja.wikipedia.org/wiki/%E3%82%A2%E3%83%AB%E3%83%95%E3%82%A1%E3%83%81%E3%83%A3%E3%83%B3%E3%83%8D%E3%83%A

Python3 & OpenCV で画像処理を学ぶ[3] 〜 トーンカーブ と LUT を理解する実装実験

1. はじめに

　本記事では、トーンカーブ（階調変換関数）を実装し、入力画像の濃淡やコントラストを補正する実験を行うことで、画素値の変換について学びます。

import numpy as np
from matplotlib import pyplot as plt
import matplotlib.gridspec as gridspec
import seaborn as sns
import cv2

# これを書くとjupyterでmatplotlibなどが出力する画像の解像度が上がる
%config InlineBackend.figure_format = 'retina'

2. トーンカーブ基礎実験

2.1 トーンカーブの基礎

　周知の通り、ディジタル画像の各画素(ピクセル)は、明暗やRGBなどの情報を表す画素値を持っています。そこで、入力の画素値 $x$ に対して出力の画素値 $y$ を対応づける関数 $f: x \rightarrow y \ (0 \leq x,y \leq 255)$ を与えることで、画像の濃淡を変換することができます。それをグラフで表したものがトーンカーブとなります。

　具体的に見ていきましょう。プロットに使用した関数のコードは長くなるので、説明の都合上この節の最後に載せます。

2.1.1 例1: $y=x$

pdr = cv2.cvtColor(cv2.imread('pandra.png',1), cv2.COLOR_BGR2RGB)
wiz = cv2.cvtColor(cv2.imread('Wizdomiot.png',1), cv2.COLOR_BGR2RGB)

def curve_1(x):
    y = x
    return y

plot_curve(curve_1,pdr)

　関数が $f(x)=x$, つまりトーンカーブのグラフが $y=x$ だと, 入力値がそのまま出力値となるので, 変化がありません。

2.1.2 例2: S字トーンカーブ

def curve_2(x):
    y = (np.sin(np.pi * (x/255 - 0.5)) + 1)/2 * 255
    return y

plot_curve(curve_2,pdr)

　このような形状のトーンカーブは俗にS字トーンカーブと言われており、 input=127 を境として、低い入力値(暗部)はより低い出力値に、高い入力値(明部)はより高い出力値に写されています。結果として暗部・明部が強調され、出力画像を見るとコントラストが高くなっています。

　ところで、コントラストの高さの定量的な特徴として、ヒストグラムが平坦であることがしばしば挙げられるようです。すなわち、コントラストが高い画像は、（写真などのように各画素値を幅広く含む画像であれば、）暗部から明部までなるべく満遍なく画素値があるということです。

　実際に、上のS字トーンカーブによる入出力のヒストグラムを見比べてみると、特に入力では少なかった暗部の方に、出力のヒストグラムの分布が広がっていることがわかります。

　このことから、トーンカーブを人手で描かなくとも、画像全体または画像を小領域に分割し、その中のヒストグラムの最小値・最大値を両端に引き延ばすようにしてコントラストを補正するアルゴリズミックな方法があります（ヒストグラム平坦化 Histogram equalization。 OpenCV-Python Tutorials のコントラスト平坦化の章を参照）。これによって、データ解析の前処理などとして、異なる照明環境で撮影された画像群のコントラストを自動で調整することができるようです。

2.1.3 例3: ガンマ変換

def curve_gamma1(x):
    gamma = 2
    y = 255*(x/255)**(1/gamma)
    return y

def curve_gamma2(x):
    gamma = 1/2
    y = 255*(x/255)**(1/gamma)
    return y

plot_curve(curve_gamma1,pdr)
plot_curve(curve_gamma2,pdr)

　このような変換は下記の関数で表され、ガンマ変換（ガンマ補正）と呼ばれています。

$y = 255\bigg(\dfrac{x}{255}\bigg)^{\frac{1}{\gamma}}$

　$\gamma > 1$ のときは1枚目（$\gamma = 2$）のように上に凸なグラフになり、すなわち全画素値が持ち上げられて画像は明るくなります。$\gamma < 1$ならばその逆で、2枚目（$\gamma = \frac{1}{2}$）のようになります。

2.1.4 使用したコード

　上記のグラフ描画のために書いた関数のコードが以下です。

def rgb_hist(rgb_img, ax, ticks=None):
    """
    rgb_img と matplotlib.axes を受け取り、
    axes にRGBヒストグラムをplotして返す
    """
    color=['r','g','b']
    for (i,col) in enumerate(color): 
        hist = cv2.calcHist([rgb_img], [i], None, [256], [0,256])
        hist = np.sqrt(hist)
        ax.plot(hist,color=col)
        
    if ticks:
        ax.set_xticks(ticks)
    ax.set_title('histogram')
    ax.set_xlim([0,256])
    
    return ax


def plot_curve(f,rgb_img):
    """
    関数 f:x->y, 画像 を受け取って
    以下のようなグラフを出す
    ----------------------------
    入力画像 | Curve | 出力画像
    histgram |       | histgram
    ----------------------------
    """
    
    fig = plt.figure(figsize=(15,5))
    gs = gridspec.GridSpec(2,3)
    x = np.arange(256)
    
    # トーンカーブを真ん中に
    sns.set_style('darkgrid')
    ax2 = fig.add_subplot(gs[:,1])
    ax2.set_title('Tone Curve')    
    ticks = [0,42,84,127,169,211,255]
    ax2.set_xlabel('Input')
    ax2.set_ylabel('Output')
    ax2.set_xticks(ticks)
    ax2.set_yticks(ticks)
    ax2.plot(x, f(x))
     
    # 入力画像を←に, 出力画像を→に
    sns.set_style('ticks')
    ax1 = fig.add_subplot(gs[0,0])
    ax1.set_title('input image →')
    ax1.imshow(rgb_img)
    ax1.set_xticks([]), ax1.set_yticks([])  # 目盛りを消す
    
    # 画素値変換。uint8で渡さないと大変なことに
    out_rgb_img = np.array([f(a).astype('uint8') for a in rgb_img]) 
    
    ax3 = fig.add_subplot(gs[0,2])
    ax3.set_title('→ output image')
    ax3.imshow(out_rgb_img)
    ax3.set_xticks([]), ax3.set_yticks([])
    
    #入力と出力のヒストグラム
    sns.set_style(style='whitegrid')
    ax4 = fig.add_subplot(gs[1,0])
    ax4 = rgb_hist(rgb_img, ax4, ticks)
    ax5 = fig.add_subplot(gs[1,2])
    ax5 = rgb_hist(out_rgb_img, ax5, ticks)
    
    plt.show()

2.2 LUT（Look Up Table）による高速化

2.2.1 LUTとは？

　先ほどまでに、入力と出力の画素値の対応づける関数を与えることで、画素値を変換しコントラストなどを補正できることを見てきました。

　ところで、一般に画像の画素数は縦×横の分だけあり、フルHD(1920×1080)であれば2073600個の画素があります。フルHDの画像の階調変換を一回行おうとすると、素朴に考えて変換の計算を2073600回行わねばならないことになります。RGB各チャンネルに変換処理を行おうとすると、その3倍になります。計算量がちょっと気がかりになりますね。

　しかしよく考えてみると、入力も出力も、取りうる値は 0~255 の整数であるため、画素値変換の対応は高々 256 通りです。とすると、いちいち変換の計算を行うのではなくて、「0~255 の各入力値に対して、出力値が 0~255 のうちどれに対応するか」を示す表をあらかじめ作っておいて、変換時はそれを参照すれば早そうです。表の作成自体は256回の計算で済みます。
　そのような表のことをLUT（Look Up Table）と呼びます。掛け算九九をいちいち計算せずに覚えておくようなものですね。以下がそのイメージです。

入力	変換	出力
0	→	0
1	→	2
…	…	…
127	→	180
…	…	…
255	→	255

2.2.2 LUTの実装

　以下に、先ほどのようにピクセルごとに計算する実装と、LUTを用いるように書き換えた実装を示します。

def naive_curve(f,rgb_img):
    # 先ほどと全く同じ変換方式
    out_rgb_img = np.array([f(a).astype('uint8') for a in rgb_img]) 
    return out_rgb_img


def LUT_curve(f,rgb_img):
    """
    Look Up Tableを LUT[input][0] = output という256行の配列として作る。
    例: LUT[0][0] = 0, LUT[127][0] = 160, LUT[255][0] = 255
    """
    LUT = np.arange(256,dtype='uint8').reshape(-1,1)
    LUT = np.array([f(a).astype('uint8') for a in LUT])
    out_rgb_img = cv2.LUT(rgb_img, LUT)
    return out_rgb_img


def simpleplot(rgb_img):
    sns.set_style('ticks')
    plt.title('output image')
    plt.imshow(rgb_img)
    plt.show()

img = naive_curve(curve_gamma1, wiz)
simpleplot(img)

　先ほどと同じ、画素ごとに変換を計算したものです。

img = LUT_curve(curve_gamma1, wiz)
simpleplot(img)

　こちらがLUTを用いた変換で、先ほどと全く同様に変換できていることがわかります。

2.2.3 LUTの速度性能

　さて、素朴に各ピクセルを変換するのと、先にLUTを作ってから変換するのとではどれだけの速度差が出るのでしょうか。計測してみましょう。

【実験1】

#先ほどのS字トーンカーブ
%timeit naive_curve(curve_2, wiz)
%timeit LUT_curve(curve_2, wiz)

52.4 ms ± 3.77 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
4.11 ms ± 46.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

【実験2】

#ガンマ変換
%timeit naive_curve(curve_gamma1, wiz)
%timeit LUT_curve(curve_gamma1, wiz)

22.5 ms ± 1.66 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)
2.95 ms ± 95.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

　明確に速度の違いが出ました。コードを見ての通り、LUTを作成するところまで含めての速度差です。LUTを用いることで、S字トーンカーブでは約12.7倍、ガンマ変換では約7.6倍もの高速化を達成することができました。

　S字トーンカーブの方の実装では三角関数を利用していたりするのでそもそもあまりよくないのですが、LUTによる変換の場合テーブル自体はすぐ作成できるため、変換に用いる関数が多少複雑なものであったりしたときも、あまり速度に差が出ないのがよいですね。

3. おわりに

3.1 今後の課題

　今回触れられなかった点として、numpy.array(dtype='uint8') は256で割った和の値であるようです。つまり、255を超えた（オーバーフローした）値は0から上がってきて、0を下回った値は255から降りてくるということです。 numpy の問題ではなく、 int 型なら一般に同じような挙動になるものと思います。
　つまり例えば、トーンカーブなどで255を上回る値を出してしまったとき、通常我々が期待するのは255でストップすることだと思われるので、どこかでそのようにクリッピングする処理を考える必要があります。階調変換関数（トーンカーブ）の段階でそのようにするのがおそらく正しそうです。

3.2 まとめと次回について

　トーンカーブとLUTの実装実験を通して、画素値を変換する際の基本的な考え方を知ることができました。

次回↓は、複数の画像の入力に対する処理について扱います。

optie.hatenablog.com

4. 参考文献

• 「ディジタル画像処理」の5章
• ヒストグラム平坦化 http://labs.eecs.tottori-u.ac.jp/sd/Member/oyamada/OpenCV/html/py_tutorials/py_imgproc/py_histograms/py_histogram_equalization/py_histogram_equalization.html#histogram-equalization

Python3 & OpenCV で画像処理を学ぶ[2] 〜 情報可視化のためのヒストグラム基礎実験

1. はじめに

　本記事では、画像の持つ特徴を定量的に知るために、その基本的な手段として、基礎的な統計量（平均, 分散など）を踏まえつつ、一次元・二次元ヒストグラムの算出・描画方法について学びます。また、HLVやL*a*b*色空間でヒストグラムを描画することで、画像データが持つ色の情報をうまく可視化することを試みます。

　前回↓

optie.hatenablog.com

はあまりコードを書きませんでしたが、今回は numpy や matplotlib で色々書いていきます。

import numpy as np
import scipy.stats as sstats
import pandas as pd
from matplotlib import pyplot as plt
import matplotlib.cm as cm
import seaborn as sns
import cv2

%matplotlib inline

2. ヒストグラム基礎実験

2.1 明度ヒストグラムと統計量の基礎

　ヒストグラムは、ある画像の画素値がどのように分布するかを示すものでした。0 から 255 までの値について(横軸)、その値を取る画素の数(縦軸)を集計してグラフにしたもの、と見ることができます。

(アルマちゃんのドット絵のヒストグラム。ドット絵なので、画素のとる値が飛び飛びなのがよくわかります)

　以下、基本的な統計量について確認しておきましょう。

　最小値[最大値]は、その名の通り、ある画像の画素値の中で最小[最大]のものです。データの分布の両端の値とも言えます。画像処理的には、最大値と最小値の差を元にコントラストが算出されたりするようです。

　平均 mean : $\mu^{2}$, 分散 variance : は、以下の式で表されます。
$k$ 個目の画素の値を$P_k$, 画像全体の画素の数を $N$ として、

平均 : $\mu = \dfrac{1}{N} \sum\limits_{k=1}^{N}P_k$

分散 :

分散の正の平方根を標準偏差 standard deviation : と言います。分散も標準偏差も、平均値から離れた画素が多いほど大きくなる統計量であるため、分布の広さ・狭さを見る指標になります。

　中央値 median は、画素値を小さい方から順に並べて、ちょうど真ん中に位置する画素値です。大きい順に並べても同じです。偶数個の場合は、真ん中にある二つの値の平均を取るようです。メディアンフィルタなどで使用されるようです。

　最頻値 mode は、画像のヒストグラムの中で一番頻度が高いものです。

　データの分布が対称な形（正規分布など）であれば、平均と分散だけでも指標として良いのですが、そうでない場合、平均は外れ値の影響を受けやすいため（例 : 平均年収はしばしば上方の少数の外れ値に引っ張られる）、関心のある目的に応じて、最頻値や中央値が分布の性質を見るために使われたりします。
　分布の偏りの方向に応じて、 最頻値≦中央値≦平均値 または 平均値≦中央値≦最頻値 という関係が成り立つようです。

　以下、ヒストグラムにこれらの統計量を表示して描画してみました。 赤い太線が平均、その両側の少し離れたところにある赤い点線が平均±標準偏差の範囲。緑の線が中央値、黄色の線が最頻値です（少々見づらいですが……）。

mean    :  114.0581962
stddev  :  65.6394397318
median  :  96.0
mode    :  70

mean    :  174.468091709
stddev  :  50.5730009287
median  :  183.0
mode    :  244

　説明の都合で後出しになりましたが、上記のヒストグラムの表示に使った関数のコードが以下です。だいたい matplotlib です。公式チュートリアルでも、ヒストグラムの表示には matplotlib を使っているようです。OpenCV の図形描画で頑張る方向性もあるようですが、普通にこちらの方が良さそうです。

def grayhist(img, stats=False):
    """
    入力 : BGR画像, 統計量の表示の有無(Optional)
    出力 : グレースケールのimshow, Histgram(+統計量)
    """
    # スタイルの設定。seabornの設定は残るため、常に最初に書いておく
    sns.set()
    sns.set_style(style='ticks')

    # プロット全体と個々のグラフのインスタンス
    fig = plt.figure(figsize=[15,4])
    ax1 = fig.add_subplot(1,2,1)
    
    sns.set_style(style='whitegrid')
    ax2 = fig.add_subplot(1,2,2)
    
    # グレースケール画像化→三重配列に戻す
    img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    img = cv2.cvtColor(img, cv2.COLOR_GRAY2RGB)
    ax1.imshow(img)
    
    img = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
    # 一次元配列化
    img = np.array(img).flatten()
    # 本来rangeは[0,255]となるはずだが、255に値が集まりすぎたので弾いた
    img = img[img!=255]

    # 軸の刻み目
    ax2.set_xticks(np.linspace(0,255,9).astype(int))
    
    # ヒストグラムを計算→色をつける
    N, bins, patches = ax2.hist(img, range=(0,255), bins=256)
    for (i,patch) in enumerate(patches):
        color = cm.gray(bins[i]/256)
        patch.set_facecolor(color)

    if stats==True: # 統計量を表示する
        mean = img.mean()
        std = np.std(img)
        median = np.median(img)
        mode = sstats.mode(img)[0][0]

        # 統計量のラインをひく
        ax2.axvline(mean, color='#d95763', linestyle='solid', linewidth=3)
        ax2.axvline(median, color='#6abe30', linestyle='solid', linewidth=2)
        ax2.axvline(mode, color='#ba9641', linestyle='solid', linewidth=2)
        ax2.axvline(mean + std, color='#d95763', linestyle='dashdot', linewidth=1)
        ax2.axvline(mean - std, color='#d95763', linestyle='dashdot', linewidth=1)

        # 統計量の説明の文字
        ax2.text(mean,N.max()*1.075, "$\mu$",color='#d95763',horizontalalignment='center')
        ax2.text(median,N.max()*1.18,"median", color='#6abe30',rotation=45)
        ax2.text(mode,N.max()*1.15,"mode",color='#ba9641',rotation=45)
        ax2.text(mean+std,N.max()*1.075, "$\mu+\sigma$", color='#d95763',horizontalalignment='center')
        ax2.text(mean-std,N.max()*1.075, "$\mu-\sigma$", color='#d95763',horizontalalignment='center')

        fig.tight_layout()
        plt.show()
        
        print("mean    : ", mean)
        print("stddev  : ", std)
        print("median  : ", median)
        print("mode    : ", mode)
        
    else:
        fig.tight_layout()
        plt.show()

実際に実行するときは以下のようになります。

# 2nd arg: 1=BGR
wiz = cv2.imread('Wizdomiot.png',1)
pdr = cv2.imread('pandra.png',1)
alma = cv2.imread('alma.png',1)


grayhist(alma)
grayhist(pdr, True)
grayhist(wiz, True)

2.2 RGB, HSV, L*a*b* 空間でのヒストグラム

　先ほどはグレースケールの明暗の値についてヒストグラムを見ました。そこで次に、色情報をヒストグラムで見ることを考えます。

　まずは基本ということで、RGBチャンネルごとのヒストグラムを見てみます。 先ほどはヒストグラムの計算にmatplotlibを用いましたが、今回はOpenCVの cv2.calcHist([img], [channel], mask_img, [binsize], ranges) を使ってみます。

wiz = cv2.cvtColor(cv2.imread('Wizdomiot.png',1), cv2.COLOR_BGR2RGB)
pdr = cv2.cvtColor(cv2.imread('pandra.png',1), cv2.COLOR_BGR2RGB)
lyrith = cv2.cvtColor(cv2.imread('Lyrith.png',1), cv2.COLOR_BGR2RGB)
sb = cv2.cvtColor(cv2.imread('SCARY BANQUET1.png',1), cv2.COLOR_BGR2RGB)
alma = cv2.cvtColor(cv2.imread('alma.png',1), cv2.COLOR_BGR2RGB)

実際に書いた関数のコードが以下です。

def rgb_hist(rgb_img):
    sns.set()
    sns.set_style(style='ticks')
    fig = plt.figure(figsize=[15,4])
    ax1 = fig.add_subplot(1,2,1)
    sns.set_style(style='whitegrid')
    ax2 = fig.add_subplot(1,2,2)
    
    ax1.imshow(rgb_img)

    color=['r','g','b']

    for (i,col) in enumerate(color): # 各チャンネルのhist
        # cv2.calcHist([img], [channel], mask_img, [binsize], ranges)
        hist = cv2.calcHist([rgb_img], [i], None, [256], [0,256])
        # グラフの形が偏りすぎるので √ をとってみる
        hist = np.sqrt(hist)
        ax2.plot(hist,color=col)
        ax2.set_xlim([0,256])

    plt.show()

実行すると、以下のようになります。

rgb_hist(wiz)

rgb_hist(pdr)

　CG制作ソフトウェアのユーザーにとっては、情報として比較的見慣れたものであると思います。
　トーンカーブなどで色調補正をするときに、いま明部にどれだけ赤成分があるか、暗部の青成分はどれくらいか、などといったことを見ながら確認したりしますね。

　次に、HSV空間とL*a*b*空間でヒストグラムを描くことを考えます。
　HSVは色相 Hue, 彩度 Saturation, 明度 Value(Brightness)の三軸からなる色空間で、L*a*b*は L* が輝度, a*が赤-緑成分, b*が黄-青成分を示す補色空間でした。(例えばa*は負値を含む軸で、プラス方向が赤、マイナス方向が緑)
　RGB空間に比べて、比較的直感的な軸であったり、人の色の認知に即した形で情報が表現されているため、より画像の色の性質をわかりやすく示してくれるかもしれません。
　以下、とりあえずHSVとL*a*b*のそれぞれのチャンネルでヒストグラムを表示する関数のコードを書いてみました。

def axes4():  
    """
    四分割したプロットを描画するときのための関数。
    axesが4つ入ったリストを返す
    """
    sns.set()
    fig = plt.figure(figsize=[15,15])
    sns.set_style(style='ticks')
    ax1 = fig.add_subplot(2,2,1)

    sns.set_style(style='whitegrid')
    ax2 = fig.add_subplot(2,2,2)
    ax3 = fig.add_subplot(2,2,3)
    ax4 = fig.add_subplot(2,2,4)

    axes = [ax1,ax2,ax3,ax4]
    return axes


def tri_hists(rgb_img, color='rgb'):
    """
    'rgb','hsv','lab'の指定を受け取って、
    「入力画像」と「各チャンネルのヒストグラム」の計4枚のプロットを描画する
    """
    axes = axes4()
    
    axes[0].set_title('image')
    axes[0].imshow(rgb_img)
    
    if color=='rgb':
        img_ch = img
        color=['R','G','B']
    elif color=='hsv':    
        img_ch = cv2.cvtColor(rgb_img, cv2.COLOR_RGB2HSV)
        color=['Hue','Saturation','Brightness']
    elif color=='lab':
        img_ch = cv2.cvtColor(rgb_img, cv2.COLOR_RGB2Lab)
        color=['L*','a*','b*']
    else:
        return

    for (i,col) in enumerate(color): # 各チャンネルのhist
        if col=='Hue':  # Hueの値域は[0,180)
            hist = cv2.calcHist([img_ch], [i], None, [180], [0,180]) 
        else:
            hist = cv2.calcHist([img_ch], [i], None, [256], [0,256])
        hist = np.sqrt(hist)
        axes[i+1].set_title(col)
        axes[i+1].plot(hist)

    plt.show()

実行結果は以下の通りです。

tri_hists(wiz,'hsv')

tri_hists(wiz,'lab')

　……いまいちよくわかりませんね。
　明るさを表すBrightnessやL*は良いとして、Hueも色相の色と対応させればわかりそうですが、「彩度が 100/255 付近に多く分布している」などと言われても解釈に困ります。

　「よくわからない」原因は色々考えられますが、いくつか挙げられるのは、彩度やa*, b*単独では最終的な「色」の想像がつきにくいことであったり、a*,b*軸は単純に馴染みが薄いことなどがあると思います。

　そこで、二軸ヒストグラムの描画を試みます。例えば、少なくとも彩度は色相とセットでないと色がわからないため、彩度軸×色相軸を組み合わせたヒストグラムを作り、グラフ上で対応する色とその頻度をセットで見てみることで、（明度を一旦無視した）「色」の頻度分布を明らかにできる可能性があります。
　この場合、 (色相, 彩度) のペアに 頻度 を対応させる二変数関数となるので、情報としては三次元となります。それをそのまま立体的な三次元グラフにしても良いのですが、隠れてしまう部分が発生するなどして見にくい可能性があるので、今回は三軸目に明度を割り当てて表現することを考えます。

　などと書いていてもイメージがつきにくいので、色相-彩度の二次元ヒストグラムを描画する関数を書いてみました。以下に示します。

def hs_hist(rgb_img):
    """
    [目的]
    「横軸に 色相 , 縦軸に 彩度 をとり、
    点(H,S)における頻度を 明度 で表現するグラフ」を作りたい。
    
    [実装]
    入力のRGB画像をHSVに変換し、(H, S)の二次元ヒストグラムを計算して(H, S, 頻度)の配列を作る。
    その配列をHSV2RGBとして変換し、RGB画像としてimshowに表示させる関数。
    """
    sns.set()
    sns.set_style(style='ticks')
    
    img_hsv= cv2.cvtColor(rgb_img, cv2.COLOR_RGB2HSV)

    # (H,S) だけにして、H,Sも分離
    hs_2d = img_hsv[:,:,0:2].reshape(-1,2)
    h = hs_2d[:,0]
    s = hs_2d[:,1]

    # ヒストグラムのbinの設定。OpenCVにおいてHの値域は0~180である（0-255に納めるためか）
    hedges = np.arange(0,180)
    sedges = np.arange(0,255)

    # 二次元ヒストグラム
    H, xedges, yedges = np.histogram2d(h,s, bins=(hedges, sedges))
    H = H.T

    # log scaleで偏りを緩和 & 正規化
    H_log = np.log(H+1)
    H_norm = H_log/H_log.max()*255

    # (H,S,頻度)の配列にするために、まずH[S]を縦[横]にリピートし、x行y列の配列にする
    x = H_norm.shape[1]
    y = H_norm.shape[0]
    hue_xy = np.repeat(xedges[:-1],y).reshape(x,y).T
    sat_xy = np.repeat(yedges[:-1],x).reshape(y,x)

    # depth方向にくっつけて、(H,S,頻度)の配列にする。uint8型でないとcvtColorが受けつけないらしい
    HS_hist = np.dstack((hue_xy, sat_xy, H_norm)).astype('uint8')
    HS_hist_im = cv2.cvtColor(HS_hist, cv2.COLOR_HSV2RGB)
    HS_hist_im = cv2.resize(HS_hist_im,(360,100))

    # 以下、元の画像をax1に、ヒストグラムをax2に表示する
    fig = plt.figure(figsize=[15,4])
    ax1 = fig.add_subplot(121)
    ax2 = fig.add_subplot(122)

    ax1.imshow(rgb_img)
    ax1.set_title('image')

    ax2.set_title('Hue-Saturation Histgram')
    ax2.set_xlabel('Hue')
    ax2.set_ylabel('Saturation')
    ax2.set_xticks(np.linspace(0,360,13))
    ax2.set_xlim(0,360)
    ax2.set_ylim(0,100)
    ax2.imshow(HS_hist_im,origin='low',interpolation='bicubic')

    plt.show() 

実行すると、以下のようになります。

hs_hist(wiz)

　縦軸に彩度、横軸に色相をとり、頻度を明度で表現した図を、左側の元画像と共に右に表示しています。こちらのWizdomiotの例でいうと、魔法陣の青〜紫のグローと、画面下部の緑色の光をヒストグラム上で見て取ることができます。

hs_hist(sb)

このSCARY BANQUETの例では、(H=270~300,S=30~50)付近に月の色の分布が目立ちます。Hue=0付近に分布しているのは、床の低彩度の赤〜橙色や本棚ですね。

hs_hist(lyrith)

　最後のLyrithの例は、元画像からは赤〜黄色系統の色が一見目立ちますが、ヒストグラムを見ると緑〜青色成分もそれなりに含まれていることがわかります。実際、画面下部の暗部などをよく見てみると、黄緑っぽい緑成分や紫がかった青成分の存在を認めることができると思います。（他人事のように解説していますが、むろん制作したのは僕なので、意図してそのようにしています）

　先ほどと比べて、「どんな色がどのように分布しているか」がわかりやすい情報が得られたのではないかと思います。左に表示した元画像と見比べても、このヒストグラムにはデータとして納得感があるのではないでしょうか。

　さて、次も同様にa*b*ヒストグラムを描いてみましょう。どうなるでしょうか。まずは書いたコードを示します。

def ab_hist(rgb_img):
    """
    基本的にhs_histと同じ
    """
    sns.set()
    sns.set_style(style='ticks')
    
    img_lab = cv2.cvtColor(rgb_img, cv2.COLOR_RGB2Lab)

    # (a,b) だけにして、a,bも分離
    ab_2d = img_lab[:,:,1:3].reshape(-1,2)
    a = ab_2d[:,0]
    b = ab_2d[:,1]

    aedges = np.arange(0,255)
    bedges = np.arange(0,255)

    H, xedges, yedges = np.histogram2d(a,b, bins=(aedges, bedges))
    H = H.T

    H_log = np.log(H+1)
    H_norm = H_log/H_log.max()*255

    # (頻度,a,b)の配列とするために、まずa[b]を縦[横]にリピートし、x行y列の配列にする。
    x = H_norm.shape[1]
    y = H_norm.shape[0]
    a_xy = np.repeat(xedges[:-1],y).reshape(x,y).T
    b_xy = np.repeat(yedges[:-1],x).reshape(y,x)

    # depth方向にくっつけて、(頻度,a,b)の配列にする。
    ab_hist = np.dstack((H_norm, a_xy, b_xy)).astype('uint8')
    ab_hist_im = cv2.cvtColor(ab_hist, cv2.COLOR_Lab2RGB)
    ab_hist_im = cv2.resize(ab_hist_im,(255,255))
    
    fig = plt.figure(figsize=[15,4])
    ax1 = fig.add_subplot(121)
    ax2 = fig.add_subplot(122)

    ax1.imshow(rgb_img)
    ax1.set_title('image')

    ax2.set_title('a*b* Histgram')
    ax2.set_xlabel('G ← a* → R')
    ax2.set_ylabel('B ← b* → Y')
    ax2.axvline(255/2, color='white', linestyle='solid', linewidth=1)
    ax2.axhline(255/2, color='white', linestyle='solid', linewidth=1)
    ax2.imshow(ab_hist_im,origin='low',interpolation='bicubic')

    plt.show()

関数を実行すると以下のようになります。

ab_hist(wiz)

　L*a*b*空間は、a*は赤-緑成分、b*は黄-青成分で、正負で補色関係になっている空間でした。すなわち、右のヒストグラムでいうと、十字の真ん中である原点(無彩色)に対して対称に分布する成分があれば、補色の成分が多く含まれていると解釈することができそうです。（軸の目盛りは無視してください……）

　上の例を見ると、基本的にa*軸のマイナス(=緑)、b*軸のマイナス(=青)成分に振れており、補色関係にあると言える成分は少ないことがわかります。わかりにくいのですが、画像上部の魔法陣の内側に赤に近い紫成分（a*軸-,b*軸+）が微妙に含まれていて、これが緑成分と補色をなしていると言えなくもないでしょうか……。

ab_hist(sb)

　こちらは、基本的にa*軸+,b*軸-の成分、すなわち紫が多いことがわかります。こちらの画像で補色を増やそうと考える場合、例えば建物内部の影になっているところの暗部に黄緑成分をトーンカーブで付加する、などといった案をヒストグラムから挙げることができそうです。

ab_hist(lyrith)

　こちらは、基本的にa*軸+,b*軸+の成分、すなわち橙色の成分がかなり多いことがわかります。彩度も幅広く分布しており、先の例に比べて色数はかなり多い方と言えるのではないでしょうか。また、a*軸+,b*軸+の橙色に対応する形で、a*軸-,b*軸-の青〜緑成分が含まれていることが確認できます（こちらは彩度など抑えめではありますが）

3. おわりに

3.1 まとめ

　今回はヒストグラムの描き方を学び、また、二次元ヒストグラムの実験を通じて、画像の特徴を可視化する方法について考えることができました。

　もちろん、「画素値の頻度を集計する」ヒストグラムは、データの可視化の手段としては比較的シンプルなものであり、「画像の特徴」を考える上で抜け落ちている重要な情報も多いと考えられます。例えば、「明度の高いピクセルが多いこと」がわかっても、それが画像中のどこに多く分布しているか？などの位置に関する情報は、ヒストグラムを通じて知ることはできません。

　しかしむしろ、上記で見てきたように、それ自体は単なる集計の可視化であるヒストグラムからでも、様々な観点からデータを見ることでわかることは想像以上に多かったというのが僕の感想です。そして今回取り上げた視点以外にも、ヒストグラムを通じて見つけられる知見はまだまだ多くありうるのではないかと思います。

3.2 今後の課題

今回至らなかった以下の点については、今後の課題としたいと思います。

1. 描画されたヒストグラムについて、値の前処理が一貫しなかったこと。

色の分布の形状を見やすくするためにと、出力を見ながら平方根をとったりlogをとったりしましたが、正当性のある補正方法が何なのかは結局わかりませんでした。

1. HSVヒストグラムの品質について。

あまり綺麗な描画ではなかったと思います。原因として考えられるのは、OpenCVがHの値を[0-180)としていること(色相の分解能が2°どまりになる)や、そもそも変換時に整数しか認められていなかったことによる情報落ちなどでしょうか。また、せっかくなので極座標に変換するなどして円形で見てみればよかったですね。

3.3 次回の方針

　次回↓は、ヒストグラムなどを用いつつ、実際に画像を操作・変換する処理の基本を学ぶ予定です。

optie.hatenablog.com

4. 参考文献

• 「ディジタル画像処理」の4章

• ヒストグラム その1: 計算して，プロットして，解析する !!! — OpenCV-Python Tutorials 1 documentation

• numpy.histogram2d — NumPy v1.12 Manual

• minus9d.hatenablog.com

Python3 & OpenCV で画像処理を学ぶ[1] 〜 色空間を工学的に理解する

　今回はあまりPython3 & OpenCVしません。色彩工学に関する基礎的な話がメインになります。

1. はじめに

　本記事は、私が「色空間についての工学的な理解を概略的に得ること」を目標とします。そのためにまず、可視光線と色認知の関係や色の三色性についての基本を復習します。次に、色の三色性に基づいて「色の見え」を定量評価する指標としてのCIE-XYZ表色系について触れ、xy色度図やL*a*b*表色系、およびこれらとRGB空間などとの関係も解説します。次に、RGB空間からHSI空間への変換について軽く述べたのちに、OpenCVでの色空間の変換を取り扱います。

　また、本記事に限りませんが、文献資料などを調べつつ自分の理解をまとめるという形を取っており、あまり厳密でない説明や曖昧な説明を含むかと思われます。もし仮に、「説明のための便宜の範疇を超えて、明らかに間違っている内容」を発見し、かつそれを指摘することが有益だと考えられる場合、ご連絡いただけると幸いです。

2. 色の基本のおさらい

　色空間の話に入る前に、まず「色」とはどのようなものだったかをおさらいしましょう。

2.1 可視光線と人が見る「色」

　通常私達が見ている「色」や「光」の物理的な正体は、電磁波 ということになります。電磁波というのは電波やX線などをも含むのですが、その中でも波長が 380nm 〜 780nm のもの（=可視光線）は人の眼で知覚することができて、そこから私たちは「色」を見ることができる、ということです。

(様々な波長を持つ電磁波のイメージ。軸が大嘘ついてますが……)

　つまり、「色」や「光」という何か特別なものが物理的世界に存在するわけではなくて、あくまで可視光線という特定の電磁波が視細胞上で化学変化を起こし、これを脳が解釈した結果として、我々は「色」を見ている……ということになります。参考までに、ニュートンがプリズムを用いた分光実験によって「光線に色は付いていない(“The rays are not coloured”)」ことを証明しているので、導出に興味のある方は調べてみてください。

　以下、「可視光線域の電磁波」を指して「光」と呼ぶことと、「色」という言葉が、「その色を知覚させる波長の電磁波」という物理的な側面と「結果として人間が見る色」という心理的な側面を併せ持つことに注意してください。

2.2. 色の三色性

可視光線スペクトル。
出展 : 可視光線 - Wikipedia

　可視光線の各波長が、それぞれの色に対応しています。
　この中で、互いに独立した3つの単色光をうまく選んで混ぜ合わせると、他の任意の色と同じ色として知覚させられることが知られています。つまり例えば、500nm付近(緑)の単色光と700nm付近(赤)の単色光を混ぜ合わせた混色光は、人間にとっては575nm付近(黄)の単色光と同じように見えるということです。緑+赤の混色光と黄の単色光は、物理的には別物であるにも関わらずです。(色の三色性)

　これは、人の視細胞上には色を弁別する受容体が三種類しかないことに起因しています。つまり、(短/中/長)波長に対して反応しやすい三種類の受容体が視細胞上に多数存在し、それぞれがどれだけ刺激されて反応したかという比率の情報によって色を判別しているため、その反応のあり方が同一になれば、 物理的には異なる光であっても、人の眼には同じ色として見えるということです。このことから、混色時に使われる三色のことを原刺激と呼びます。
　関連して、いわゆる太陽や蛍光灯などの白く見える光とは、(誤解を恐れずに言えば、)これらの三種類の受容体がバランスよく反応するような、複数の波長の光を含む可視光線だと考えることができます。

太陽光に含まれる波長ごとの強さを示す分光エネルギー分布。可視光線領域が最も高くなっている。
出展 : 太陽光 - Wikipedia

　任意の色が三色だけの混合で表せるという性質が色の三色性ですが、これを利用した代表的なものが、ディスプレイのRed, Green, Blueや、印刷のCyan, Magenta, Yellow(+K)ですね。RGBは闇にR,G,Bの光を足していったもの(=加法混色)で、CMYは光にC,M,Yフィルタを通したもの、すなわち白色光からR,G,Bを引き去ったものになります(=減法混色)。「R,G,Bを引き去る」とは、例えば白色光をYellowフィルタに通すと黄色く見えますが、これは白色光がフィルターを通るときにBlue成分が吸収されて無くなった結果そう見える、ということです。

RGBの加法混色と、CMYの減法混色。CMYKのKとはKey Plateの略で、C,M,Yのみでは完全な黒にならないために追加されるインク。
出展 : 色の３原色 (T. Fujiwara)

3. XYZ表色系 - 色の見えを規定する

　光は電磁波であって色は認知上の産物であること、三色の混色で任意の色を表現できることがわかりました。以上の事実を踏まえて、「色」を定量的な指標で表すことを考えましょう。以下、色を表示するための体系である表色系について触れます。

3.1 CIE-XYZ表色系

　さて、先ほど「任意の単色は、三色（=三刺激）の混合によって等色させられる」ことを述べました。この事実から、実験によって「ある色を等色するのに必要な三刺激値」を求めることができそうです（等色実験）。実際そのようにして、「ある波長 λ の光に対して、等色に必要な三刺激の値を返す三組の関数」である等色関数を構成することができます。

CIE-XYZ等色関数。
出展 : CIE 1931 色空間 - Wikipedia

　そうして、国際照明委員会（CIE）により、可視光線380 ~ 780nmで定義された等色関数の組が CIE-XYZ表色系 です。実験では、三刺激（RGB）の値は部分的に負の値をとる*1ようなのですが、そのような関数は実用上問題があるので、全て正になるように調整されています。また、CIE-XYZ表色系では、波長ごとの人間の「明るさ」の知覚(比視感度曲線)も$\overline{y}(\lambda)$で近似できるように設定されており、$Y$だけで色の明るさを規定できます。*2

　実際の光は特定波長のみを持つ単色光というわけではなく、しばしば複数の波長を持つ合成光（太陽光など）であるわけですが、どのような色光であっても、その光が含む波長を示す分光エネルギー分布 $L(\lambda)$ が与えられれば、以下のような積分で三刺激値を求めることができます。*3

$ X = \int_{380}^{780} \overline{x}(\lambda)L(\lambda) d\lambda $

$ Y = \int_{380}^{780} \overline{y}(\lambda)L(\lambda) d\lambda $

$ Z = \int_{380}^{780} \overline{z}(\lambda)L(\lambda) d\lambda $

3.2 xy色度図

　次に、こうして得られる三刺激値 $X,Y,Z$ の相対比 $x,y,z$ を考えます。すなわち、 $X,Y,Z$ それぞれの刺激値は、三刺激値全体 $X+Y+Z$ の中でどれくらいの割合を占めるか？というようなことを考えます。（例えば、 刺激値 $X $ の大きさが全体に占める割合が80%なら、$(x,y,z)=(0.8, 0.1, 0,1)$ のような形になります）そうすると、刺激値の相対比 $x,y,z$ は、以下のような式で表されることがわかると思います。

$ x = \dfrac{X}{X+Y+Z} $

$ y = \dfrac{Y}{X+Y+Z} $

$ z = \dfrac{Z}{X+Y+Z} $

また、これらの式とその意味（ $x,y,z$ が相対比であること）から、

$ x + y + z = 1 $

であることもわかるかと思います。つまり、$ z=1-x-y $ なので、$x,y$ の値が決まれば $z$ もただ一つの値に決まるということです。このことから $(x,y)$ を二次元座標で表現すると、「明るさを無視した、三刺激値の比率と色の対応」の図を作ることができ、それが以下のxy色度図となります。$(x,y)=(1,0),(0,1)$の二点を結ぶ対角線の右上側は、$x+y$だけですでに1を超える領域のため存在しないことに注意してください。($x,y,z$は非負なので)

xy色度図。
出展 : CIE 1931 色空間 - Wikipedia

　ここには、人間の知覚できる(明るさの違いを除いた)全ての色が、歪んだ楕円を切断したような領域の中に表されています。このような形状になる理由は、例えば$(x,y)=(1.0,0,0)$となるような光は存在しない、すなわち $X$ の刺激値だけを上げるような光は物理的に存在せず、光が$X$の刺激値を上げるならば必ず隣接する $Y$ にも影響するからであり、したがって、この形状は物理的に取りうる $(x,y)$ の領域であるということができます。
　ディスプレイに表示する色の観点では、LEDの三原色R,G,Bの発色の割合を $ (r,g,b)=(1.0,0.0,0.0),(0.0,1.0,0.0),(0.0,0.0,1.0) $ のようにしたとき、その三点を三刺激値 $(X,Y,Z)$ およびその比率 $(x,y)$ に対応させてこのxy色度図にプロットし、それを結んだ三角形がディスプレイに表示できる色の範囲ということになります。

3.3 CIE-L*a*b*表色系

　前節で、等色実験に基づいて光の波長と三刺激値の対応を与えるCIE-XYZ表色系について見ましたが、$(X,Y,Z)$ の値をそのまま三次元座標として捉えた空間は、等色関数の重なりなどの影響から、色空間としては歪んでおり、そのままでは色に対して何か操作を考えるには困難がありそうです。
　そこで、$X,Y,Z$ に適当な変換を施して、我々の「空間」の直感に合う直交座標の三次元色空間を再現する目的で設計されたのが、以下のような CIE-L*a*b表色系 です。

L*a*b*表色系の図。
出展 : Lab 表色系 (T. Fujiwara)

　 $(X,Y,Z) \rightarrow (L^*, a^*, b^*) $ とするための具体的な変換は、以下のようなものです。

$ L^* = 116 \left( \dfrac{Y}{Y_0} \right)^{ \frac{1}{3} } - 16 $

$ a^* = 500 \left( \left( \dfrac{X}{X_0} \right)^{\frac{1}{3}} - \left( \dfrac{Y}{Y_0} \right)^{\frac{1}{3}} \right) $

$ b^* = 200 \left( \left( \dfrac{Y}{Y_0} \right)^{\frac{1}{3}} - \left( \dfrac{Z}{Z_0} \right)^{\frac{1}{3}} \right) $

　ここで $ X_0, Y_0, Z_0 $ は、対象と同一の照明条件の下の標準白色に対する三刺激値です。要するに、照明によって色の見え方が違うため、その照明における「白色」の三刺激値をもとに基準化をしているということです。
　式はやや複雑ですが、$a^*$ はCIE-XYZ図のグラフ的に $X$ と重なる $Y$ の影響を引いて排除したもの、すなわち（負値を含む）赤成分で、$b^*$ も同様に黄成分だと考えればよいと思います。$L^*$ は明るさです。

　CIE-L*a*b*表色系の具体的なメリットとしては、二色 $ (L^*_1, a^*_1, b^*_1 ), ( L^*_2, a^*_2, b^*_2 ) $ の色の差 $\Delta$ を、以下のように三次元空間における普通の距離(ユークリッド距離)として簡単に計算できることです。

$ \Delta = \sqrt{ \left( L^*_1 - L^*_2 \right)^2 + \left( a^*_1 - a^*_2 \right)^2 + \left( b^*_1 - b^*_2 \right)^2 } $

3.4 XYZ表色系とRGB値との関係

　さて、我々が馴染み深いのは（XYZなどではなく）RGB画像データであるわけですが、いわゆる $(R,G,B)$ と先程までのXYZ表色系などとはどのような関係にあるのでしょうか。
　察しがついているかもしれませんが、 $(R,G,B)$ の値それ自体は具体的な色を表しているわけではなくて、ディスプレイなどによる出力を経て初めて具体的な色になります。したがって、モニターがデータとして $(R,G,B)$ の値を受け取ってどのように発光するか（実際の色の見え=XYZ表色系 にどう対応づけられる発光になるのか）などといったことによって、データとしてのRGB値が表す具体的な色は異なってきます。
　 そこで、 $(R,G,B) \rightarrow (X,Y,Z)$ という具体的な対応づけを規格化したものが、sRGBやAdobe RGBなどだということになります。これは、前述のxy色度図において三角形で表されるものです。モニターの発色による $(R,G,B) \rightarrow (X,Y,Z)$ の対応と、sRGBなどの規格が定める $(R,G,B) \rightarrow (X,Y,Z)$ の対応とどれだけ近いかが、「正しい発色」を考える上で問題になっていることです。

4. RGBとHSI

　RGBや、色相,彩度,明度からなるHSI(HLSとも)の概念は周知であると思います。 R,G,Bの三軸による表現はあまり人間にとって直観的ではないので、代わりに色相,彩度,明度というわかりやすい属性で色を扱おうという話ですね。
　しかしながら、その相互変換の方法はいくらか存在するようです。六角錐モデル、双六角錐モデル、円柱モデルなどがありますが、今回はその中でも、人間の直感により近いという意味で精密とされる双六角錐モデルについて以下で解説します。

4.1. 双六角錐モデルによるRGB-HSI変換

はじめに、双六角錐モデルの模式図を以下に示します。次に、数式的な解説を行うので、見比べながら確認してください。

双六角錐モデルの模式図。白/黒に近い色は彩度が低いという点が特徴である。他のモデルでは、白と原色の赤などが同じ彩度1.0とされることもある。
出展 : 色モデル (T. Fujiwara)

　まず、明度 $I$ は以下のようになります。

$ I_{max} = max\{R,G,B\} $

$ I_{min} = min\{R,G,B\} $

$ I = \dfrac{I_{max} + I_{min}}{2} $

こちらは大丈夫でしょうか。$(R,G,B)$ の中で一番大きい値と一番小さい値の平均が明度であると言っていますね。

　次に、彩度 $S$ です。

$ S = \left\{ \begin{array}{ll} \dfrac{I_{max} - I_{min}}{I_{max} + I_{min}} & (I \leqq 0.5) \\ \dfrac{I_{max} - I_{min}}{2-(I_{max} + I_{min})} & (I < 0.5) \end{array} \right. $

$R=G=B$ のときの色がグレースケール、すなわち $I_{max} = I_{min}$ となり彩度 $S=0$ となること。グレーを境に、明度が上がって白に近づくほど、また明度が下がって黒に近づくほど、分母が大きくなり彩度は小さくなること。また、彩度が最大になるのは、$I_{max} = 1$ かつ $I_{min} = 0$ , すなわち $I=0.5$ なるときであることを確認してください。

　最後に、色相 $H$ です。 $H$ は角度であることに注意して、

$r = \dfrac{I_{max}-R}{I_{max} - I_{min}}$

$g = \dfrac{I_{max}-G}{I_{max} - I_{min}}$

$b = \dfrac{I_{max}-B}{I_{max} - I_{min}}$

とおくと、これらの文字を使って、以下のように表されます。

$ H = \left\{ \begin{array}{ll} \dfrac{\pi}{3}(b-g) & (R=I_{max}のとき) \\ \dfrac{\pi}{3}(2+r-b) & (G=I_{max}のとき) \\ \dfrac{\pi}{3}(4+g-r) & (B=I_{max}のとき) \end{array} \right. $

　例えば、 $R=I_{max}$ であるときは色相として赤成分が一番強いということに注意すると、係数の $\dfrac{\pi}{3}$ および定数項と場合分けの条件から、円周をR,G,Bを中心に三等分して、それぞれの支配領域の中で残り二色の比率を考えているものと見ることができそうです。
$I_{max} = I_{min}$ (グレースケール)であるとき、分母がゼロとなるので色相は定義されません。

5. OpenCVによる色変換実験

　ようやくOpenCVの話です。が、OpenCVでは cvtColor(image, flag) 関数の第二引数に「何から何への変換か」を示すflag (例: cv2.COLOR_RGB2GRAY ) を与えるだけで変換がされたimageが返ってきます。
　以下、もぺもぺを用いて実験してみましょう。各色空間での値をRGBの値であると言い張ってmatplotlibに表示させ、ピクセル配列の一部と共に表示しました。

from matplotlib import pyplot as plt
import numpy as np
import cv2

%matplotlib inline


IMAGE_PATH = '../0.IntroToOpenCV/mopemope_title.jpg'

img_gbr = cv2.imread(IMAGE_PATH, 1)
img_rgb = cv2.cvtColor(img_gbr, cv2.COLOR_BGR2RGB)

img_gray = cv2.cvtColor(img_rgb, cv2.COLOR_RGB2GRAY)
img_xyz = cv2.cvtColor(img_rgb, cv2.COLOR_RGB2XYZ)
img_lab = cv2.cvtColor(img_rgb, cv2.COLOR_RGB2LAB)
img_hls = cv2.cvtColor(img_rgb, cv2.COLOR_RGB2HLS)

spaces = [img_rgb, img_gray, img_xyz, img_lab, img_hls]
names = ['RGB', 'Grayscale', 'XYZ', 'L*a*b*', 'HLS']

for space, name in zip(spaces, names):
    plt.title('mopemope' + name)
    if name=='Grayscale':
        # GRAYの場合(width, height)の二次元配列となるため、三次元配列に変換しておく
        space = cv2.cvtColor(space, cv2.COLOR_GRAY2RGB)
    plt.imshow(space)
    plt.show()
    print('----- part of pixels -----')
    print(space[0])
    print('====================')

    ----- part of pixels -----
    [[217 255 216]
     [217 255 216]
     [217 255 216]
     ..., 
     [255 255 255]
     [255 255 255]
     [255 255 255]]
    ====================

---

    ----- part of pixels -----
    [[239 239 239]
     [239 239 239]
     [239 239 239]
     ..., 
     [255 255 255]
     [255 255 255]
     [255 255 255]]
    ====================

---

    ----- part of pixels -----
    [[220 244 240]
     [220 244 240]
     [220 244 240]
     ..., 
     [242 255 255]
     [242 255 255]
     [242 255 255]]
    ====================

---

    ----- part of pixels -----
    [[246 109 143]
     [246 109 143]
     [246 109 143]
     ..., 
     [255 128 128]
     [255 128 128]
     [255 128 128]]
    ====================

---

    ----- part of pixels -----
    [[ 59 236 255]
     [ 59 236 255]
     [ 59 236 255]
     ..., 
     [  0 255   0]
     [  0 255   0]
     [  0 255   0]]
    ====================

また、変換に使うflagですが、flags = [i for i in dir(cv2) if i.startswith('COLOR_')] のようにして探すことができます。 以下では、デジタル完結の作業にはあまり関係がないと思われるベイヤーフィルターを弾いて、冒頭6文字 'COLOR_' を除いたものを変換元ごとに表示するコードです。

import re

flags = [i for i in dir(cv2) if i.startswith('COLOR_')]
flags = [flag for flag in flags if not re.match(r'.*(Bayer|BAYER).*',flag)]


A=''

for flag in flags:
    if A == re.split(r'\w(2)\D',flag)[0]:
        print(flag[6:], end=', ')
    else:
        A = re.split(r'\w(2)\D',flag)[0]
        print('\n')
        print(flag[6:], end=', ')

以下、その出力です。

    BGR2BGR555, BGR2BGR565, BGR2BGRA, BGR2GRAY, BGR2HLS, BGR2HLS_FULL, BGR2HSV, BGR2HSV_FULL, BGR2LAB, BGR2LUV, BGR2Lab, BGR2Luv, BGR2RGB, BGR2RGBA, BGR2XYZ, BGR2YCR_CB, BGR2YCrCb, BGR2YUV, BGR2YUV_I420, BGR2YUV_IYUV, BGR2YUV_YV12, 
    
    BGR5552BGR, BGR5552BGRA, BGR5552GRAY, BGR5552RGB, BGR5552RGBA, 
    
    BGR5652BGR, BGR5652BGRA, BGR5652GRAY, BGR5652RGB, BGR5652RGBA, 
    
    BGRA2BGR, BGRA2BGR555, BGRA2BGR565, BGRA2GRAY, BGRA2RGB, BGRA2RGBA, BGRA2YUV_I420, BGRA2YUV_IYUV, BGRA2YUV_YV12, 
    
    COLORCVT_MAX, 
    
    GRAY2BGR, GRAY2BGR555, GRAY2BGR565, GRAY2BGRA, GRAY2RGB, GRAY2RGBA, 
    
    HLS2BGR, HLS2BGR_FULL, HLS2RGB, HLS2RGB_FULL, 
    
    HSV2BGR, HSV2BGR_FULL, HSV2RGB, HSV2RGB_FULL, 
    
    LAB2BGR, LAB2LBGR, LAB2LRGB, LAB2RGB, 
    
    LBGR2LAB, LBGR2LUV, LBGR2Lab, LBGR2Luv, 
    
    LRGB2LAB, LRGB2LUV, LRGB2Lab, LRGB2Luv, 
    
    LUV2BGR, LUV2LBGR, LUV2LRGB, LUV2RGB, 
    
    Lab2BGR, Lab2LBGR, Lab2LRGB, Lab2RGB, 
    
    Luv2BGR, Luv2LBGR, Luv2LRGB, Luv2RGB, 
    
    M_RGBA2RGBA, 
    
    RGB2BGR, RGB2BGR555, RGB2BGR565, RGB2BGRA, RGB2GRAY, RGB2HLS, RGB2HLS_FULL, RGB2HSV, RGB2HSV_FULL, RGB2LAB, RGB2LUV, RGB2Lab, RGB2Luv, RGB2RGBA, RGB2XYZ, RGB2YCR_CB, RGB2YCrCb, RGB2YUV, RGB2YUV_I420, RGB2YUV_IYUV, RGB2YUV_YV12, 
    
    RGBA2BGR, RGBA2BGR555, RGBA2BGR565, RGBA2BGRA, RGBA2GRAY, RGBA2M_RGBA, RGBA2RGB, RGBA2YUV_I420, RGBA2YUV_IYUV, RGBA2YUV_YV12, RGBA2mRGBA, 
    
    XYZ2BGR, XYZ2RGB, 
    
    YCR_CB2BGR, YCR_CB2RGB, 
    
    YCrCb2BGR, YCrCb2RGB, 
    
    YUV2BGR, YUV2BGRA_I420, YUV2BGRA_IYUV, YUV2BGRA_NV12, YUV2BGRA_NV21, YUV2BGRA_UYNV, YUV2BGRA_UYVY, YUV2BGRA_Y422, YUV2BGRA_YUNV, YUV2BGRA_YUY2, YUV2BGRA_YUYV, YUV2BGRA_YV12, YUV2BGRA_YVYU, YUV2BGR_I420, YUV2BGR_IYUV, YUV2BGR_NV12, YUV2BGR_NV21, YUV2BGR_UYNV, YUV2BGR_UYVY, YUV2BGR_Y422, YUV2BGR_YUNV, YUV2BGR_YUY2, YUV2BGR_YUYV, YUV2BGR_YV12, YUV2BGR_YVYU, YUV2GRAY_420, YUV2GRAY_I420, YUV2GRAY_IYUV, YUV2GRAY_NV12, YUV2GRAY_NV21, YUV2GRAY_UYNV, YUV2GRAY_UYVY, YUV2GRAY_Y422, YUV2GRAY_YUNV, YUV2GRAY_YUY2, YUV2GRAY_YUYV, YUV2GRAY_YV12, YUV2GRAY_YVYU, YUV2RGB, YUV2RGBA_I420, YUV2RGBA_IYUV, YUV2RGBA_NV12, YUV2RGBA_NV21, YUV2RGBA_UYNV, YUV2RGBA_UYVY, YUV2RGBA_Y422, YUV2RGBA_YUNV, YUV2RGBA_YUY2, YUV2RGBA_YUYV, YUV2RGBA_YV12, YUV2RGBA_YVYU, YUV2RGB_I420, YUV2RGB_IYUV, YUV2RGB_NV12, YUV2RGB_NV21, YUV2RGB_UYNV, YUV2RGB_UYVY, YUV2RGB_Y422, YUV2RGB_YUNV, YUV2RGB_YUY2, YUV2RGB_YUYV, YUV2RGB_YV12, YUV2RGB_YVYU, 
    
    YUV420P2BGR, YUV420P2BGRA, YUV420P2GRAY, YUV420P2RGB, YUV420P2RGBA, 
    
    YUV420SP2BGR, YUV420SP2BGRA, YUV420SP2GRAY, YUV420SP2RGB, YUV420SP2RGBA, 
    
    YUV420p2BGR, YUV420p2BGRA, YUV420p2GRAY, YUV420p2RGB, YUV420p2RGBA, 
    
    YUV420sp2BGR, YUV420sp2BGRA, YUV420sp2GRAY, YUV420sp2RGB, YUV420sp2RGBA, 
    
    mRGBA2RGBA, 

6. おわりに

　本記事では、可視光線から復習し、等色実験から色の見えを定量的に規定するXYZ表色系やxy色度図、それを元にしたL*a*b*色空間や、RGB&HSIの変換と工学的な立ち位置を概観することができたと思います。また、変換の数式的な理解をざっくり得られたことで、各色空間の間の関係も掴みやすくなったのではないかと思います。
　次回↓は、具体的な画像を用いて、画像の統計量に関して取り扱おうと思います。

optie.hatenablog.com

また、一応考えている今後の予定としては、以下の通りです。

[done] 0. 導入
[done] 1. 色空間 
2. 画像の統計量
3. 画素ごとの濃淡変換
4. 空間フィルタリング
5. 周波数フィルタリング
6. テクスチャ・領域処理
7. パターンと図形検出
8. パターン認識

　省略・割り込み・入れ替えなど起こるかもしれませんが、概念的には「画素ベースの操作」→「画素の位置関係ベースの操作」→「画素の周波数ベースの操作」→「特徴量抽出」→「パターン認識」という流れで学習を進めることがが念頭にあります。今回もそうなのですが、理論的な教科書としては暫定的に「ディジタル画像処理」（CG-ARTS協会）に従いつつ、OpenCVによる実験を行いながら学んでいく感じになると思います。

7. 参考文献

• 「ディジタル画像処理」の3章

• CIE 1931 色空間 - Wikipedia

• 第1回 大事なのは"正しい色"を表示できること――液晶ディスプレイの「色域」を理解しよう | EIZO株式会社

• http://w3.kcua.ac.jp/~fujiwara/infosci/lab.html

*1:一般に光を足し合わせると明るくなる=彩度が落ちるので、高彩度の Cyan などは Red をマイナスしないと物理的には等色できない、といったことがあるようです。

*2:というよりもともと、波長ごとに「明るさ」の感じ方を示す比視感度曲線は、緑成分=中波長付近で最大になるような形をしており、中波長付近に反応するM錐体細胞の感度曲線とほぼ一致しているようです。これは、物理的に同じ強さの光であれば、緑色を最も明るく感じることを示唆しています。グラフ上の三刺激の中でも $\overline{y}(\lambda)$ の最大値が最も低いのも、緑色を最も強く感じるために、刺激値としては少ない量で済むことを表していると考えることもできそうです。

*3:補足 : 「入力光の中の ある波長 $\lambda$ nmの強さに対して、その刺激値を求める」という計算を、入力光に含まれる 380 ~ 780nm の範囲の波長全てに対して行い、それらを足し合わせるというイメージです。結果として、入力光に対する刺激値を得ます

Python3 & OpenCV で画像処理を学ぶ[0] 〜 OpenCVの導入

　これから何回かに分けて、Python3でOpenCVライブラリを使用しながら画像処理を僕が学んでいきたいと思います。 目標としては、「OpenCVでできること&その書き方を概略的に把握し、未知の実装も適当に調べれば書ける状態になること」といちおう定義しておきます。

1. 概要

1.1 動機

　画像処理ができたら……嬉しい！
　画像処理についてきちんと学んだ経験はありませんが、今のところ「特徴量抽出、(画素|領域|周波数)単位での(操作|変換|フィルタリング)、パターン認識」といったあたりに関心を持っています。

1.2 なぜPythonか

　Pythonが僕にとって一番慣れた言語であること、numpyなどとの連携ができること、jupyter が使えるとうれしい等といった理由があります。

1.3 これからやること

　今回は、OpenCVの初歩的な使い方を手を動かしながら覚えます。 具体的には、公式チュートリアルの Gui Features in OpenCV を(適宜省略などしつつ)進めます。
　適宜コメントをつけて理解を深めますが、コードはほぼ写経です。

　また、Python, numpy, matplotlibの書き方に関する基本的な知識と、CG映像ソフトのユーザー程度のディジタル画像に関する知識(R,G,Bチャンネルの数値を持つピクセルという概念や、ディジタル画像はピクセルの配列だと理解していることなど)は仮定します。

今回使用した環境は以下の通りです。

- MacOS High Sierra 10.13.3
- Python 3.6.2
- OpenCV 3.3.0_3
- Jupyterlab 0.27.0

諸々のインストールが完了しているところから話を始めます。

　それでは、やっていきます。必要ライブラリをインポートし、サンプルとして使用する素材のパスを変数にしておきます。

import cv2
import numpy as np
from matplotlib import pyplot as plt

%matplotlib inline  #jupyterなので

IMAGE_PATH = 'mopemope_title.jpg'
MOV_PATH = 'Lyrith_170906.mp4'

2. 画像の入出力と表示

　画像は cv2.imread(path, mode) 関数で読み込めるようです。画像なので、ピクセルの二重配列が返り値です。 カラーの場合は、横 * 縦 * [B,G,R]の三重配列になります。OpenCVでは、 Red, Green Blue ではなく Blue, Green, Red の順番のようです。 可視光線の短波長(紫外線方向)から長波長(赤外線方向)の順の並びなので、こちらの方が自然ということでしょうか。

画像の書き込みはcv2.imwrite(path, image) となります。

　「画像を読み込んでopenCVのウィンドウで表示し、sキーが入力されれば画像を保存する」コードが以下になります。

# 第二引数は1:RGB, 0:GrayScale, -1:RGBA
img = cv2.imread(IMAGE_PATH, 0)

# 第二引数 cv2.WINDOW_NORMAL : ウィンドウをリサイズできるようにする
cv2.namedWindow('image', cv2.WINDOW_NORMAL)

cv2.imshow('image',img)
k = cv2.waitKey(0) # any key inputを無限に待つ

if k==27:  #ESC
    cv2.destroyAllWindows()
elif k == ord('s'):  #sキー
    cv2.imwrite(OUT_IMAGE_PATH,img)
    cv2.destroyAllWindows()

出力:

無事に、グレースケールで読み込んだもぺもぺのタイトル画像が表示されました。

次に、画像をmatplotlibで表示してみます。

img = cv2.imread(IMAGE_PATH, 0)
plt.imshow(img, cmap='gray',interpolation='bicubic')
# plt.xticks([]), plt.yticks([])  #目盛りをなくす
plt.show()

出力:

　問題なく表示できました。 このあたりは直感的でわかりやすいですね。

3. 映像の入出力と表示

　映像の入力の場合、 VideoCapture オブジェクトを呼び、.read() メソッドでフレームを一枚づつ読んでいくようです。 .read() メソッドは返り値として、フレームと「フレームがあるかどうか」の真偽値(終了判定) を持つため、 これを利用してループを回し、フレームを一枚づつ読み書きするという形になります。

　映像を出力する際は、 VideoWriter オブジェクトを使用します。インスタンスを作るときに、cv2.VideoWriter(filename, codec, fps, frame_size) という形で書き出し形式を引数で与えます。

　コーデックを選択するときなのですが、 Video Codecs by FOURCC - fourcc.org にある4バイトの識別コードを cv2.VideoWriter_fourcc() 関数に渡した返り値を用いるようです。今回 ANIM を利用しましたが、先ほどのページを見ても説明がなく、ANIMが実際何なのかはよくわからないのですが、これを.movに与えたものがQuickTime Animationに相当するのではないかと推察されます(雑)。また、cv2.VideoWriter_fourcc(*'ANIM') ないし cv2.VideoWriter_fourcc('A','N','I','M')と書く必要があるそうです。*'ANIM' という表現は何なのでしょうか、ポインタではないような……？

　VideoWriterオブジェクトの .write(frame) メソッドに、実際に保存したいフレームを次々に与えることで書き出されるようです。

　以下は、リリスの動画を読み込んで上下反転させて表示し、指定したファイルに書き込むコードです。

cap = cv2.VideoCapture(MOV_PATH)

#codec  http://www.fourcc.org/codecs.php
fourcc = cv2.VideoWriter_fourcc(*'ANIM')

#VideoWriterObj. args: name, codec, fps, frame size
out = cv2.VideoWriter('output.mov',fourcc, 20.0, (1920,1080)) 

while(cap.isOpened()):
    #ret : True/False(EOF判定), frame : 各フレーム
    ret, frame = cap.read()
    if ret:  #終了フレームまで
        frame = cv2.flip(frame,0)  #上下反転
        
        out.write(frame)

        cv2.imshow('frame', frame)
        
        #waitKey(ms)を挟み,　キー入力 q がなければ引き続きwhileループを継続
        if cv2.waitKey(1) & 0xFF ==ord('q'):
            break
    else:
        break

cap.release()
out.release()
cv2.destroyAllWindows()

出力:

上下が反転した迷宮リリスが表示されました。

4. 基本的な描画機能

　図形の描画です。はじめに背景となる画像(配列)を用意し、そこに図形が表現されるように配列の数値を書き換えていくという流れになります。 試しに黒平面に対角線を引いてみます。

# 0 埋めされたwidth * height * (B,G,R)の三重配列 (黒平面)
img = np.zeros((512,512,3), np.uint8)

# 線を引く args: array, start, end, (B,G,R), thickness 
# 返り値はimg配列の値を操作したような配列
img = cv2.line(img, (0,0), (511,511), (255,0,0), 5)

　これをmatplotlibで表示してみます。 先ほども述べたように、OpenCV上ではあくまで(B,G,R)として扱われているので、(R,G,B)で解釈するmatplotlibのためにcv2.cvtColor(img, cv2.COLOR_BGR2RGB)と変換して渡してあげます。

#matplotlibは(R,G,B)で解釈するので、変換して渡す
plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
plt.show()

出力:

このように描画されました。

以下、基本的な図形の描画です。長方形、円、楕円はコードのコメントを見てください。

cv2.polylines(img, [pts], isClosed, color, thickness)は与えられた座標を結ぶ線分を描画するのですが、 座標の与え方として、 のnumpy配列を、さらに[]でくくって一次元リストの要素にしたものでないといけないようです。これは何故なのでしょう。

# 長方形 args: array, ↖︎, ↘︎, (B,G,R), thickness 
img = cv2.rectangle(img, (350, 350), (400, 500), (0,155,155), 10)

# 円 args: array, 中心, 半径, (B, G,R), thickness. 
# 閉じた図形のthicknessに-1 を与えるとfillとなる
img = cv2.circle(img, (447, 63), 63, (0,0,155), -1)

#楕円 args: array, 中心, (長半径, 短半径), angle, startAngle, endAngle, Color, thickness
# θ ∈ [startAngle, endAngle] の弧を作り、その後にangleで図形自体を回転させる
img = cv2.ellipse(img, (256,256), (100, 50), 45, 30, 210, (255,155,0), 10)

pts = np.array([[10,5], [200,30], [70, 200], [50, 100]], np.int32)
pts = pts.reshape((-1,1,2))
#折れ線 args: array, 点群, isClosed, Color, thickness
img = cv2.polylines(img, [pts], False, (0,255,255), 5)

テキストは以下のように書けます。HERSHEYというフォントのファミリーしか使えない？

font = cv2.FONT_HERSHEY_SCRIPT_COMPLEX

# args : array, text, ↙︎, fontFace, font_Scale, Color, thickness, linetype(アンチエイリアス)
img = cv2.putText(img, 'OptieCV', (10,300), font, 4, (255,255,255), 2, cv2.LINE_AA)

以上の描画結果を表示します。

plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
plt.show()

出力:

先ほどの対角線の上に、長方形,円,楕円,折れ線,テキストが描画されました。

こうして、OpenCVで基本的な図形を描画する方法をわかることができました。

5. GUIトラックバーの作成

　ユーザーがインタラクティブに値を操作できるスライダーを実装します。 まず cv2.namedWindow(window_name) 関数でウィンドウを定義し(名前をつけて)、cv2.createTrackbar(Trackbar_label, window_name, 初期値, 最大値, callback_func) 関数でウィンドウにトラックバーを割り当てます。callback_func はバーを動かすたびに呼ばれる関数です。何もないときは pass するだけの関数を渡しておけば良いです。

以下、R,G,Bそれぞれの値を操作可能なスライダーを与えて、その色を表示するウィンドウを作るコードです。

def nothing(x):
    pass

img = np.zeros((300, 512, 3), np.uint8)

cv2.namedWindow('image')

# args: label, window_name, min, max, callback_func
cv2.createTrackbar('R', 'image', 0, 255, nothing)
cv2.createTrackbar('G', 'image', 0, 255, nothing)
cv2.createTrackbar('B', 'image', 0, 255, nothing)

switch = '0 : OFF \n 1 : ON'
cv2.createTrackbar(switch, 'image', 0, 1, nothing)

while(1):
    cv2.imshow('image',img)
    k = cv2.waitKey(1) & 0xFF
    if k==27:
        break
    
    r = cv2.getTrackbarPos('R','image')
    g = cv2.getTrackbarPos('G','image')
    b = cv2.getTrackbarPos('B','image')
    s = cv2.getTrackbarPos(switch,'image')
    
    if s==0:
        img[:] = 0
    else:
        img[:] = [b, g, r]

cv2.destroyAllWindows()

出力:

　……チュートリアル通りのコードのはずなのですが、目盛りがおかしいことと、switchの改行が見切れていること、バーの現在の値が表示されていませんね。 一応、ちゃんと動いてはいるので今回は良しとします。

　動作も重いので、インタラクティブに操作をしたい場合は jupyter 側で ipywidgets などを導入した方が良いかもしれません。

6. まとめ

　今回は、OpenCVによる画像/映像の入出力と取り扱い、図形の描画に関して導入することができました。 今のところOpenCVのGUI操作を利用する予定はないので、次回以降の方針として、画像への操作・変換・フィルタリングを重点的に見ていきたいと思います。

7. 参考文献

Gui Features in OpenCV — OpenCV-Python Tutorials 1 documentation

https://docs.opencv.org/2.4/modules/core/doc/drawing_functions.html#ellipse

[Python]PythonでOpenCVを使う (基本編) - Qiita

このブログは？

はじめに

　Optieです。ブログを始めました。

　本記事は、最初の投稿として、「このブログが一体何であるのか」をまとめることとします。 そのために、「なぜ書くのか」「何を書くのか」という点について簡潔に記します。

なぜ書くのか

　ブログ開設の主な理由として、「アウトプットの練習」があります。

　元々、思考の言語化に対しては少し苦手意識があります。特に、人に伝える前提で整った文章を書くことは、僕は得意としていません。 その大きな原因の一つとして、アウトプット量が少ないことが考えられます。そこで、本ブログを通じて、言語的なアウトプットの練習をできればと考えました。 要するに経験値稼ぎなのですが、いま考えていることや新しく知ったこと等を、自力で言語化していく作業を行なっていきたいです。

何を書くのか

　まだ具体的な記事などは決めていないのですが、「考えたことと、その実装」「新しく知ったこと」「活動の記録」といったことを書ければと考えています。 例えば、作品の制作時に考えたことや、技術的な思いつきや練習などで手を動かしたこと、読んだ本や大学で学んだこと、今月やったことなど。 アウトプット自体が目的でもあるので、映像や制作の話に限らず、関心のあることを広く書いていけたらと思っています。

おわりに

　正直これ継続できるかわからないのですが、できれば週一くらいで書けたらいいなあとか。まあこういう感じなので、温かく見守って頂ければと……。