1. おことわり

この記事は筆者の備忘録として書いており、一般の方が読むことを想定していません。未来の自分が読み返して理解できるように、できるだけ分かりやすく書くつもりですが、前提となる説明が省略されている箇所が多数あります。ご了承下さい。

2. 背景と目的

Jzazbz Color Space を使い BT.709 色域の Gamut Boundary を求めようとしたが、特定の Hue では筆者の使用していたアルゴリズム(Method A)で破綻が発生した
前回の記事では破綻を防ぐ暫定的なアルゴリズム（Method B）を使用したが、計算コストが高く困っていた
そこで今回は新しいアルゴリズム (Method C) を考えることにした

3. 結論

計算時間の削減を実現した新しいアルゴリズム (Method C) を考案した
従来の方法 (Method B) と比較して計算時間が200～1000倍程度早くなった。
- Method B の発想が小学生レベルだったのが原因だが…

4. 理論

4.1. Gamut Boundary の簡単な説明

Gamut Boundary とは Color Space の外枠である。CIELAB および Jzazbz Color Space での BT.709 色域の Gamut Boundary の例を動画1, 動画2 に示す。


動画1. CIELAB の Gamut Boundary	動画2. Jzazbz の Gamut Boundary

Gamut Boundary の算出は筆者が趣味で行っている色彩工学の勉強では重要である。例えば筆者はこれまで以下の用途で Gamut Boundary を使用してきた。

① Gamut Mapping の検討 (過去記事へのリンク)
② 目視で表示デバイスの品質を判断するテストパターンの作成 (過去記事へのリンク)

その一方で、高精度に Gamut Boundary を求めることは意外と難しいことが分かってきた。以後で、Gamut Bundary の扱い方、および算出アルゴリズムについて簡単に説明していく。

4.2. Gamut Boundary の扱い方（2D LUTとして扱う方法の紹介）

Gamut Boundary の計算アルゴリズムの話に入る前に、Gamut Boundary のデータ形式や扱い方を簡単に紹介する（なお、これらは筆者が独自に考えたものであり一般的な方法かどうかは不明である）。

まず Gamut Boundary のデータ形式について説明する。筆者は Gamut Boundary を L* (Lightness)、 C* (Chroma)、h* (Hue angle) の3次元データとして扱っている。また、Gamut Boundary は 2D LUT として保存している。 2D LUT を使う理由は Gamut Boundary の計算は（アルゴリズムにも依るが）時間がかかるためである。事前に計算して 2D LUT 化すると簡単な線形補間計算だけで任意の Gamut Boundary を得られるため、とても便利である。

2D LUT の入力と出力は以下の通りである。

入力： L* (Lightness) と h* (Hue angle)
出力：C* (Chroma) を含む L*C*h*

具体例を示そう。以下はサンプルとして作成した 8x8 の 2D LUT である。サンプルコードにあるように lut[lightness_idx][hue_idx] (※1) の入力に対して [ 71.43 74.89 102.86] の出力 (※2)を得ることができる。

※1 lightness_idx=5 は L*=71.43 を、hue_idx=2 は 102.86° を意味する
※2 (L*, C*, h*) = (71.43, 74.89, 102.86) である

>>> lut = np.load("./lut/2dlut_8x8.npy")
>>> print(lut)
[[[   0.      0.      0.  ]  # L*, C*, h*
  [   0.      0.     51.43]
  [   0.      0.    102.86]
  [   0.      0.    154.29]
  [   0.      0.    205.71]
  [   0.      0.    257.14]
  [   0.      0.    308.57]
  [   0.      0.    360.  ]]

 [[  14.29   36.24    0.  ]
  [  14.29   27.58   51.43]
  [  14.29   21.39  102.86]
  [  14.29   23.25  154.29]
  [  14.29   13.68  205.71]
  [  14.29   17.6   257.14]
  [  14.29   75.05  308.57]
  [  14.29   36.24  360.  ]]

 [[  28.57   53.34    0.  ]
  [  28.57   51.21   51.43]
  [  28.57   38.1   102.86]
  [  28.57   34.22  154.29]
  [  28.57   20.14  205.71]
  [  28.57   25.91  257.14]
  [  28.57  110.45  308.57]
  [  28.57   53.34  360.  ]]

 [[  42.86   70.43    0.  ]
  [  42.86   68.79   51.43]
  [  42.86   50.42  102.86]
  [  42.86   45.19  154.29]
  [  42.86   26.59  205.71]
  [  42.86   34.21  257.14]
  [  42.86  115.27  308.57]
  [  42.86   70.43  360.  ]]

 [[  57.14   79.96    0.  ]
  [  57.14   85.49   51.43]
  [  57.14   62.65  102.86]
  [  57.14   56.16  154.29]
  [  57.14   33.05  205.71]
  [  57.14   42.52  257.14]
  [  57.14   85.2   308.57]
  [  57.14   79.96  360.  ]]

 [[  71.43   48.07    0.  ]
  [  71.43   58.16   51.43]
  [  71.43   74.89  102.86]
  [  71.43   67.13  154.29]
  [  71.43   39.5   205.71]
  [  71.43   45.37  257.14]
  [  71.43   55.89  308.57]
  [  71.43   48.07  360.  ]]

 [[  85.71   21.82    0.  ]
  [  85.71   24.03   51.43]
  [  85.71   87.12  102.86]
  [  85.71   78.1   154.29]
  [  85.71   45.95  205.71]
  [  85.71   22.36  257.14]
  [  85.71   27.47  308.57]
  [  85.71   21.82  360.  ]]

 [[ 100.      0.      0.  ]
  [ 100.      0.     51.43]
  [ 100.      0.    102.86]
  [ 100.      0.    154.29]
  [ 100.      0.    205.71]
  [ 100.      0.    257.14]
  [ 100.      0.    308.57]
  [ 100.      0.    360.  ]]]

>>> lightness_idx = 5
>>> hue_idx = 2
>>> print(lut[lightness_idx][hue_idx])
[  71.43   74.89  102.86]

それでは（やや無謀だが）この 2D LUTを使い L*=70 の a-b 平面をプロットしてみる。やり方は簡単である。lightness_idx = 5 時のデータをプロットするだけである。

プロットした結果を図2に示す。

f:id:takuver4:20210903195752p:plain:w500 — 図2. 8x8 の 2D LUT を使って Gamut Boundary をプロットした様子

当然のことながら、LUT の Grid数が少なすぎて正しいプロットが行えていない。では Grid数を 16x16, 32x32, 64x64 と増やしていこう。結果を図3～図5 に示す。


図3. 16x16 の 2D LUT を使用	図4. 16x16 の 2D LUT を使用	図5. 64x64 の 2D LUT を使用

Grid数を増やしていくと理想的な形に近付くことが分かる。筆者は図1 をプロットする際は 1024x1024 の 2D LUT を使用した（※3）。

※1 厳密には LUT値をそのまま使うのでは無く線形補間して使用した。本記事では線形補間の詳細は省略する（単純に線形補間しているだけなので）。

さて、ここまでは特定の L*値に対する2次元の Gamut Boundary のプロットをしてきた。次は L*値も変化させて3次元の Gamut Boundary をプロットする。 8x8, 16x16, 32x32, 64x64 でプロットした結果を以下に示す。


動画3. 8x8, 16x16, 32x32, 64x64 でプロットした様子

4.3. アルゴリズム

続いて、Gamut Boundary の求め方を紹介する。ここでは Method A, Method B, Method C の3種類を紹介する。

4.3.1. Method A

（ここから先の文章は不自然な改行が多くなります。ご了承下さい。ブログの数式表示の仕様です。）

Method A は反復処理で Gamut Boundary を探索する方法である。概要を図6 に示す。

f:id:takuver4:20210912150516p:plain:w600 — 図6. Method A で Gamut Boundary を探索する様子 (k=0)

まず、Gamut Boundary より十分外側に点 $\displaystyle P_{j, k=0}$ を置く。

ここで j は 2D LUT の hue angle の index (※4) を、 k は反復処理の index を意味する。また、 $\displaystyle P_{j, k=0}$ の Chroma を

$\displaystyle C^{*}_{j, 0}$ とする。

この条件で点 $\displaystyle P_{j, k+1}$ は次のようにして求める。

もしも、点 $\displaystyle P_{j, k}$ が Gamut Boundary の外側だった場合は

$\displaystyle \frac{C^{*}_{j, 0}}{2^{k+1}}$ だけ Chroma を減少させる。

もしも、点 $\displaystyle P_{j, k}$ が Gamut Boundary の内側だった場合は

$\displaystyle \frac{C^{*}_{j, 0}}{2^{k+1}}$ だけ Chroma を増加させる。

点 $\displaystyle P_{j, k=0}$ の場所にも依るが、24回程度の反復処理で十分な精度の Gamut Boundary を求めることができる。

※4 図6 は 8x8 の 2D LUT に対して i=1 でプロットした。hue angle 値は 360 / (8 - 1) = 51.4° となる。

4.3.2. Method B

Method A が上手く機能するのは hue angle の Gamut Boundary が 1つしか存在しない 場合のみである。実は CIELAB color space や Jzazbz color space は L*, C*, h* の極座標系で表現しようとすると、Gamut Boundary が1つ以上存在する場合がある。CIELAB の例を図7 に、Jzazbz の例を図8, 図9 に示す。

説明	図
図7 CIELAB の例
図8 Jzazbz の例（全体）
図9 Jzazbz の例（拡大）

こういった状況の場合、Method A では Gamut Boundary が異常な形状になることがある。

筆者はこの状況について以下の方針で対応することにした。

完璧な Gamut Boundary の算出は断念する
- L*, C*, h* の 2D-LUT を使う場合、原理的に完璧な表現が不可能
- L*, C*, h* の 2D-LUT を使わない方法は個人的にやりたくなかった
色相のズレを最小化するため、Gamut Boundary が複数個見つかった場合は Chroma が最小のものを選択する

この方針に従って考えた方法が Method B である。

Method B は以下の記事に詳細がある。本記事では一部の文言だけ引用して細かい説明は省略する。

trev16.hatenablog.com

1. 任意の Hue の Chroma-Lightness 平面を Chroma 方向に N_c個、
    Lightness方向に N_l個に分割する（N_c=20, N_l=10 の例を図5に示す）
2. 各ブロックの LCH値を RGB値に変換し、0 ≦ R ≦ 1 かつ 0 ≦ G ≦ 1 かつ 0 ≦ B ≦ 1 を
    満たすブロックを有効ブロックとする（有効ブロックを黄色にした例を図6に示す）
3. 各 Lightness値に対して Chromaが大きくなる方向にブロックを走査し、
    最初に無効ブロックとなる直前のブロックを、その Lightness値の Gamut Boundary とする
   （Gamut Boundary を黒色にした例を下図に示す）
4. これを Hue を少しずつ変更しながら 0～360°の範囲で行う

f:id:takuver4:20210923091833p:plain — 図10. Method B の説明の図

Method B によって問題の解決が可能となった。しかし新たに別の問題が生じた。計算コストが非常に高い、という問題である。筆者が求める精度で計算をしようとすると 2D LUT の生成に 100時間以上かかることが分かった。筆者は様々な色域、様々な輝度値に対する 2D LUT を作成したかったのだが、Method B を使った場合は現実的な計算時間では終わらないことが判明した。

そこで考えたのが、本記事の主題である Method C である。

4.3.3. Method C

Method C は簡単に言うと Method A と Method B を組み合わせた方法である。最初に Method B で Gamut Boundary を探索する領域を決めて、その後に Method A で精度の高い Gamut Bundary を求める 2-pass 方式である。

CIELAB での例を図11、図12 に示す。


図11 Method C (全体)	図12 Method C ( $\displaystyle Q_{j, n=4}$ 付近を拡大)

図11 では Method B を使用して粗い精度で Gamut Boundary を求めている（この例では $\displaystyle N=8$ ）。ここでは赤色で塗った箇所が Gamut Boundary と判定されている。

また、このとき $\displaystyle Q_{j, n=N-1}$ の Chroma を $\displaystyle C^{\ast}_{Q, j}$ , これを N個の点で分割した Chroma を $\displaystyle C^{*}_{P, j}$ とする。

$\displaystyle C^{\ast}_{P, j} = \frac{C^{\ast}_{Q, j}}{N-1}$ である。

図12 では Method A を使用して高い精度で Gamut Boundary を求めている。まず Method B を使用して求めた点を $\displaystyle P_{j, k=0}$ とする。

この条件で点 $\displaystyle P_{j, k+1}$ は次のようにして求める。

もしも、点 $\displaystyle P_{j, k}$ が Gamut Boundary の外側だった場合は

$\displaystyle \frac{ C^{\ast}_{P, j} }{2^{k}}$ だけ Chroma を減少させる。

もしも、点 $\displaystyle P_{j, k}$ が Gamut Boundary の内側だった場合は

$\displaystyle \frac{ C^{\ast}_{P, j} }{2^{k}}$ だけ Chroma を増加させる。

Method B のサンプル数 $\displaystyle N$ にも依るが、20回程度の反復処理で十分な精度の Gamut Boundary を求めることができる。

4.4. 結果の確認

今回、Method C を使って CIELAB Color Space, Jzazbz Color Space の Gamut Boundary の 2D LUT を作成した。作成した 2D LUT は以下の方法で確認を行った。

4.4.1. 目視による確認

まずは大きな間違いが無いかを目視で確認した。作成した 2D LUT を使用して CIELAB, Jzazbz の双方で a-b(az-bz)平面, C*-L*(CzJz)平面をプロットして明らかな破綻が無いことを確認した。確認用に作成した動画を以下に示す。


目視確認用に作成した動画

確認していたところ、想定通りの破綻と想定外の破綻があった。代表的なものを図13 と図14 に示す。


図13 a-b平面での破綻例	図14 CzJz 平面での破綻例

図13 は想定通りの破綻である。前述の通り CIELAB Color Space は極座標系での完璧な表現が不可能である。今回のアルゴリズムでは Gamut Boundary は Chroma が小さくなるよう設計したので、設計通りの結果になっている。

図14 は想定外の破綻であった。詳細はここでは記載しないが、幾つかの追加デバッグを行った結果、筆者はこの破綻を 2D LUT の HUE 方向のサンプル数不足が原因だと結論づけた。今回作成した 2D LUT は、hue angle 方向のサンプル数を 4096 にしたのだが、250°～260° は変化が激しく線形補間して使用すると破綻が生じてしまう。サンプル数 16384 まで増やすことも考えたが、以下の理由により対応しないことにした。