Color Transfer between Imagesを実装した

2020年5月24日追記
タイトルを「Color Transfer between Imagesを実装したが論文のような結果にはならなかった」から「Color Transfer between Imagesを実装した」に変更しました．色がおかしくなっていた原因はnp.uint8を行う前に0~255の範囲に値を限定していなかったからです．

本文

機械学習を利用しない古典的な画像の色変換の論文Color Transfer between ImagesをPyTorchで実装しました．この論文は入力画像（source image）の色の雰囲気を参照画像（target image）に近づける画像処理アルゴリズムです．
PyTorchを使った理由はtensorの扱いに慣れたかったからです．

実行結果

入力画像	参照画像	自分の結果	論文中の結果	追加の結果

自分の結果は入力画像の色を参照画像に近づけた結果です．追加の結果は入力画像と参照画像を逆にしたものです．
イラストに対してもやってみました．

入力画像	参照画像	結果1	結果2

入力画像を参照画像に近づけるたのが結果1，逆にしたものが結果2です．

具体的な画像処理アルゴリズム

入力画像の色の雰囲気を参照画像に近づけるのはRGB色空間ではなく，$l \alpha \beta$色空間という色空間でします．まずRGB色空間から$l \alpha \beta$色空間への計算式について書き，その後入力画像を参照画像の色に近づける計算式を書きます．

sRGBからlαβ色空間への変換

論文ではRGB色空間からXYZ色空間への変換式を利用していますが，現在ではsRGBを利用しているようなので，sRGBからXYZ色空間への変換式を利用しました．sRGBからXYZ色空間への変換式は教科書を参照しました．
$[0, 255]$の値を持つ8bit sRGBを$255$で割る事で$[0, 1]$の範囲に変換します．R,G,Bの各チャネルの全ての値に対して次の変換をします．

f_{1}\left(c\right)=

\left\{ 
    \begin{array}{ll}
\frac{c}{12,92} & \text { for } c \leq 0.04045 \\
\left(\frac{c+0.055}{1.055}\right)^{2.4} & \text { for } c>0.04045
\end{array}

\right.

その後，各ピクセル値に対して次の行列計算して，XYZ色空間へと変換します．

\left[
\begin{array}{l}
X \\ Y \\ Z
\end{array}
\right]
=M_{\mathrm{RGB}}^{-1} \left[\begin{array}{l}R \\ G \\ B\end{array}\right]

ただし

M_{\mathrm{RGB}}^{-1}=
\left[
\begin{array}{ccc}
0.412453 & 0.357580 & 0.180423 \\ 
0.212671 & 0.715160 & 0.072169 \\ 
0.019334 & 0.119193 & 0.950227
\end{array}
\right]

です．
XYZ色空間の値をLMS色空間の値に変換します．

\left[\begin{array}{c}L \\ M \\ S\end{array}\right]
=\left[\begin{array}{ccc}
0.3897 & 0.6890 & -0.0787 \\ 
-0.2298 & 1.1834 & 0.0464 \\ 
0.0000 & 0.0000 & 1.0000
\end{array}\right]
\left[\begin{array}{l}X \\ Y \\ Z\end{array}\right]

LMS色空間ではデータが歪んでいる（skew）ためlogをとります．

\begin{array}{c}
\mathbf{L}&=\log L\\
\mathbf{M}&=\log M\\
\mathbf{S}&=\log S
\end{array}

次の計算を行うことで$l \alpha \beta$色空間の値に変換します．

\left[\begin{array}{c}l \\ \alpha \\ \beta\end{array}\right]=
\left[\begin{array}{ccc}
\frac{1}{\sqrt{3}} & 0 & 0 \\ 
0 & \frac{1}{\sqrt{6}} & 0 \\ 
0 & 0 & \frac{1}{\sqrt{2}}
\end{array}\right]
\left[\begin{array}{ccc}
1 & 1 & 1 \\
1 & 1 & -2 \\
1 & -1 & 0
\end{array}\right]
\left[
\begin{array}{c}\mathbf{L} \\ \mathbf{M} \\ \mathbf{S}\end{array}
\right]

一連の計算はクラスConversion_Srgb_LMSのsrgb_to_labに実装しています．

入力画像の色の雰囲気を参照画像に近づける

入力画像の各チャネルの平均を引きます．

\begin{aligned}
&l^{*}=l-\langle l\rangle\\
&\alpha^{*}=\alpha-\langle\alpha\rangle\\
&\beta^{*}=\beta-\langle\beta\rangle
\end{aligned}

次に入力画像の分散で割り，参照画像の分散を掛けることで参照画像の色のばらつきを与えます．

\begin{aligned}
&l^{\prime}=\frac{\sigma_{t}^{l}}{\sigma_{s}^{l}} l^{*}\\
&\alpha^{\prime}=\frac{\sigma_{t}^{\alpha}}{\sigma_{s}^{\alpha}} \alpha^{*}\\
&\beta^{\prime}=\frac{\sigma_{t}^{\beta}}{\sigma_{s}^{\beta}} \beta^{*}
\end{aligned}

最後に，$l^{\prime},\alpha^{\prime},\beta^{\prime}$に参照画像の平均を加えます．$l\alpha\beta$色空間内で新しい画像を生成することができました．
この計算は関数color_transfer内に実装してあります．

lαβ色空間からsRGBへの変換

$l \alpha \beta$空間内の画像を以下の計算でsRGB空間内に戻します．
この計算はクラスConversion_Srgb_LMS内のlab_to_srgbに実装してあります．

\left[\begin{array}{c}
\mathbf{L} \\
\mathbf{M} \\
\mathbf{S}
\end{array}\right]=\left[\begin{array}{ccccc}
1 & 1 & 1 \\
1 & 1 & -1 \\
1 & -2 & 0
\end{array}\right]\left[\begin{array}{ccc}
\frac{\sqrt{3}}{3} & 0 & 0 \\
0 & \frac{\sqrt{6}}{6} & 0 \\
0 & 0 & \frac{\sqrt{2}}{2}
\end{array}\right]\left[\begin{array}{c}
l \\
\alpha \\
\beta
\end{array}\right]

\begin{array}{c}
L &=  \mathrm{e}^{\mathbf{L}}\\
M &= \mathrm{e}^{\mathbf{M}}\\
S &= \mathrm{e}^{\mathbf{S}}
\end{array}

\left[\begin{array}{l}R \\ G \\ B\end{array}\right]
=M_{\mathrm{RGB}}
\left[\begin{array}{ccc}
0.3897 & 0.6890 & -0.0787 \\ 
-0.2298 & 1.1834 & 0.0464 \\ 
0.0000 & 0.0000 & 1.0000
\end{array}\right]^{-1}
\left[\begin{array}{c}L \\ M \\ S\end{array}\right]

M_{\mathrm{RGB}}=\left[\begin{array}{rrr}
3.240479  & -1.537150 & -0.498535 \\
-0.969256 & 1.875992  & 0.041556 \\
0.055648  &-0.204043  & 1.057311
\end{array}\right]

最後にRGB各チャネルの全ての値に以下のガンマ補正をすることでsRGB画像が生成できます．

f_{1}(c)=\left\{\begin{array}{ll}
12.92 \cdot c & \text { for } c \leq 0.0031308 \\
1.055 \cdot c^{1 / 2.4}-0.055 & \text { for } c>0.0031308
\end{array}\right.

ソースコード

import torch
from PIL import Image
import numpy as np
from pathlib import Path


class Conversion_Srgb_LMS:
    def __init__(self):
        srgb_to_xyz = torch.tensor([
            [0.412453, 0.357580, 0.180423],
            [0.212671, 0.715160, 0.072169],
            [0.019334, 0.119193, 0.950227]])

        xyz_to_lms = torch.tensor([
            [0.3897, 0.6890, -0.0787],
            [-0.2298, 1.1834, 0.0464],
            [0.0000, 0.0000, 1.0000]])
        self.srgb_to_lms = torch.mm(xyz_to_lms, srgb_to_xyz)
        self.lms_to_srgb = self.srgb_to_lms.inverse()
        m1 = torch.tensor(
            [[1, 1, 1], [1, 1, -2], [1, -1, 0]], dtype=torch.float32)
        vec = torch.rsqrt(torch.tensor([3.0, 6.0, 2.0]))
        m2 = torch.diag(vec)
        self.LMS_to_lab_matrix = torch.einsum("mc, cd ->md", m2, m1)

        m1 = m1.T
        vec = torch.div(torch.sqrt(torch.tensor(
            [3.0, 6.0, 2.0])), torch.tensor([3, 6, 2], dtype=torch.float32))
        m2 = torch.diag(vec)
        self.lab_to_LMS_matrix = torch.einsum("mc, cd->md", m1, m2)

    def srgb_to_lab(self, srgb_img):
        """
        input srgb_img tensor [0, 1]
        return lab_img tensor
        """
        srgb_img = torch.where(srgb_img <= 0.04045, srgb_img/12.92,
                               torch.pow((srgb_img+0.055)/1.055, 2.4))
        lms_img = torch.einsum("whc, mc -> whm", srgb_img, self.srgb_to_lms)
        LMS_img = torch.log(lms_img)
        lab_img = torch.einsum("whc, mc-> whm", LMS_img,
                               self.LMS_to_lab_matrix)
        return lab_img

    def lab_to_srgb(self, lab_img):
        """
        input lab tensor
        return srgb tensor[0, 1]
        """
        LMS_img = torch.einsum("whc, mc -> whm", lab_img,
                               self.lab_to_LMS_matrix)
        lms_img = torch.exp(LMS_img)
        srgb_img = torch.einsum("whc, mc -> whm", lms_img, self.lms_to_srgb)
        srgb_img = torch.where(srgb_img <= 0.0031308, 12.92 *
                               srgb_img, 1.055 * torch.pow(srgb_img, 1/2.4) - 0.055)
        srgb_img = torch.clamp(srgb_img, min=0, max=1)

        return srgb_img


def color_transfer(source_img, target_img):
    """
    source_img: lab image tensor
    target_img: lab image tensor

    return:参照画像（target img）の色に変換した入力画像(source img)
    """
    source_mean = torch.mean(source_img, dim=[0, 1])
    target_mean = torch.mean(target_img, dim=[0, 1])
    source_std = torch.std(source_img, dim=[0, 1], unbiased=False)
    target_std = torch.std(target_img, dim=[0, 1], unbiased=False)
    new_lab_img = torch.div(target_std, source_std) * \
        (source_img - source_mean) + target_mean
    return new_lab_img



def main(source_path, target_path, save_dir):
    source_img = np.asarray(Image.open(source_path)) / 255
    # αチャネルがあった場合は削除
    if source_img.shape[2] == 4:
        source_img = source_img[:, :, :3]
    source_img = torch.from_numpy((source_img.astype(np.float32)))
    target_img = np.asarray(Image.open(target_path)) / 255
    if target_img.shape[2] == 4:
        target_img = target_img[:, :, :3]
    target_img = torch.from_numpy(target_img.astype(np.float32))
    conversion = Conversion_Srgb_LMS()
    source_lab_img = conversion.srgb_to_lab(source_img)
    target_lab_img = conversion.srgb_to_lab(target_img)
    new_lab_img = color_transfer(source_lab_img, target_lab_img)
    new_srgb_img = conversion.lab_to_srgb(new_lab_img)
    img = Image.fromarray(np.uint8(255*new_srgb_img.cpu().float().numpy()))
    save_name = Path(save_dir).joinpath(
        "source=" + Path(source_path).stem + "target=" + Path(target_path).stem + ".png")
    img.save(save_name)

if __name__ == "__main__":
    import sys
    args = sys.argv
    if len(args) != 4:
        print("color_tranfer.py 入力画像へのパス 参照画像へのパス 保存するディレクトリへのパス")
        print("と実行してください")
        sys.exit()
    main(args[1], args[2], args[3])

論文中のRGBからXYZ色空間へ変換行列の導出

論文中のRGBからXYZ色空間へ変換行列は次の通りです．

\left[\begin{array}{l}
X \\
Y \\
Z
\end{array}\right]=
\left[\begin{array}{lll}
0.5141 & 0.3239 & 0.1604 \\
0.2651 & 0.6702 & 0.0641 \\
0.0241 & 0.1228 & 0.8444
\end{array}\right]\left[\begin{array}{l}
R \\
G \\
B
\end{array}\right]

この行列を導出します．
国際通信連合（International Telecommunication Union, ITU）がCIE 1931RGBからXYZ色空間への変換する規格CIE XYZitu601-1 (D65)を決めたようです．決められた行列は次の通りで，この値を修正して上記の行列を導出します．

M_{\mathrm{itu}}=\left[\begin{array}{ccc}
0.4306 & 0.3415 & 0.1784 \\
0.2220 & 0.7067 & 0.0713 \\
0.0202 & 0.1295 & 0.9394
\end{array}\right]

表記を簡単にするために

M_{\mathrm{itu}}=
\left[\begin{array}{ccc}
a&b&c\\
d&e&f\\
g&h&d
\end{array}\right]

と置きます．
ここで，白色($R=G=B=1$)で$x = \frac{X}{X+Y+Z}=0.333\cdots,y=\frac{Y}{X+Y+Z}=0.333\cdots$となるようにしたいです．
しかし，そのまま$M_{\mathrm{itu}}$を利用した場合，$R=G=B=1$のとき

\left[
\begin{array}{l}
X \\ Y \\ Z
\end{array}
\right]
=
M_{\mathrm{itu}}
\left[\begin{array}{c}
1 \\
1 \\
1
\end{array}\right]
=
\left[\begin{array}{c}
0.9505\\
1.0000\\
1.0891
\end{array}\right]

となりとなり，$x =\frac{X}{X+Y+Z} \neq 0.333\cdots,y=\frac{Y}{X+Y+Z} \neq 0.333\cdots$と先ほどの条件を満たすことができません．

そこで

M_{\mathrm{itu}}x=
\left[\begin{array}{c}
1 \\
1 \\
1
\end{array}\right]

となる$x$を求めると

x = 
\left[\begin{array}{c}
1.1940\\
0.9483\\
0.9081
\end{array}\right]

となります．
$x^{T}$を$M_{\mathrm{itu}}$の各行にかけて

\left[\begin{array}{ccc}
1.1940a & 0.9483b& 0.9081c\\
1.1940d & 0.9483e& 0.9081f\\
1.1940g & 0.9483h& 0.9081d
\end{array}\right]
=
\left[\begin{array}{lll}
0.5141 & 0.3239 & 0.1604 \\
0.2651 & 0.6702 & 0.0641 \\
0.0241 & 0.1228 & 0.8444
\end{array}\right]

となり論文中の行列を求めることができ，$R=G=B=1$とき，$X=Y=Z=1$となっています．

$M_{\mathrm{itu}}は$色彩工学に乗っている下記のCIE 1931 RGB色空間からXYZ色空間への変換行列と異なるため，どのようにして求められたか気になりました．また，白色($R=G=B=1$)で$x = \frac{X}{X+Y+Z}=0.333\cdots,y=\frac{Y}{X+Y+Z}=0.333\cdots$となるようにしたい理由も完全に理解できませんでした．

\left[
\begin{array}{c}
X \\
Y \\
Z
\end{array}\right]
=\left[
\begin{array}{ccc}
2.76883 & 1.75171 & 1.13014 \\
1.0 & 4.59061 & 0.06007 \\
0.0 & 0.05651 & 5.59417
\end{array}\right]
\left[\begin{array}{l}
R \\
G \\
B
\end{array}\right]

lαβ色空間を利用する理由

$l \alpha \beta$色空間は人間の視覚モデル，ヘーリングの反対色説をモデル化した色空間であるため，画像処理に適している．この空間のデータはlogスケールのため値の分散が小さい．また，対照的でもあるらしい．各軸の相関が弱いため，各軸で異なる処理を行ってもアーティファクトが起こりにくい．

試してもダメだったこと

sRGBからXYZ色空間への変換式ではなく，論文中の値をそのまま使っても，同じように色がおかしくなることが置きました．また，sRGBの変換行列を$R=G=B=1$のとき，$X=Y=Z=1$となるよう行列を修正しても生成画像の色がおかしくなりました．この修正時には簡単にするためガンマ補正をしませんでした．
上手くいかなかった原因はnp.uint8を行う前に0~255の範囲に値を限定していなかったからです．

感想と疑問点

色って難しい．画像変換系のコードを書くとき，どうやってデバックすれば良いのだろう？
コードが間違っているのか，それとも入力画像のペアが悪かったのかどうすれば気づけるのだろうか．
スクリーンショットで画像を作成したからおかしな結果になった可能性はあるのだろうか．
LMS色空間ではデータが歪んでいる（skew）ためlogをとるという，意味もよく分からなかった．
$l \alpha \beta$色空間への変換式の行列のかっこの大きさが揃っていなくて気持ち悪い．綺麗な書き方ってどうやるんですかね．
$l \alpha \beta$色空間でデータが対照的になるのは何故だろうか．