1. はじめに

この記事は「Report ITU-R BT.2407 の Annex2 を実装」シリーズの3回目です。かなり間が空きましたが失踪せずに少しずつ実装を進めていました。無事に実装が完了したので結果を報告したいと思います。

2. 目的

BT.2407 Annex2 に記載の BT.2020 to BT.709 の Gamut Mapping を実装する。

自身の環境を踏まえて以下の点を意識した実装とする。

①Limited Range の信号前提ではなく Full Range の信号を前提とした処理とする。

BT.2407 Annex2 は NHKが考案した方式ということもあり、現在のTV放送に使用されている Limited Range の信号（940Lv～1019Lv の Overwhite 領域も含む）で扱えるように、CIELAB色空間を拡張して処理を行っている。

一方で自分は基本的に RGB444 Full Range の信号処理にしか興味が無いため、Limited Range を考慮した処理は行わないこととした。

②DCI-P3 to BT.709 の変換も視野に入れた実装とする

BT.2407 Annex2 は BT.2020 to BT.709 変換に最適化したアルゴリズムとなっている。が、自分は BT.2020 以外にも DCI-P3 の信号も扱うことが多い。そのため将来的には本実装を DCI-P3 to BT.709 変換にも応用したいと考えている。

こうした背景もあり BT.2407 Annex2 の Hue Mapping は今回の実装には盛り込まないこととした。理由は BT.2020 to BT.709 に特化したマジックナンバーが複数あり、DCI-P3 to BT.709 変換への応用が困難だと判断したからである（※1）。

※1偉そうに書いていますが、別の言い方をすると自分の頭ではマジックナンバーの意味を正しく理解することが無理だったということです。おバカでごめんなさい…。

3. 結論

どうにか実装が完了した。加えて簡単に扱えるように 3DLUT化も行った。結果を図1～図6に示す。

図1は効果確認用のテストパターンである。水平方向が Hue の変化、垂直方向が Chroma の変化を意味する。また Chroma が BT.709 の Gamut Boundary を超えたブロックは左上に黒いポチを付けてある。本アルゴリズムは BT.709 の外側の色にのみ変換を行う方法のため、内側の色に影響が出ていないことの確認をしやすいよう黒ポチを付けた。

図2～図6 は図1を BT.2020 to BT.709 変換した結果である。より詳しい説明は以下を参照。（はてなブログの画像圧縮を嫌って Google Drive から画像を引っ張ってきています。表示に少し時間がかかる場合があります）

図1：BT.2020空間での評価用テストパターン（※2）
図2：BT.2020空間での Gamut Mapping の実行結果（※2）
図3：BT.2020 to BT.709 の Gamut Mapping したものを BT.709 色域で表示したもの（※3 ）
図4：筆者にとっての従来手法である 3x3 の Matrix で変換したものを BT.709 色域で表示したもの。図3との比較用（※3 ）
図5：図1と図2をアニメーションで切り替え（※2）
図6：図3と図4をアニメーションで切り替え。筆者にとっての従来手法と BT.2470 Annex2 との比較（※3 ）

※2 BT.2020カバー率100%の表示デバイスでご確認下さい
※3 一般的な BT.709色域の表示デバイスでご確認下さい



図1. 効果確認用テストパターン(BT.2020色域, Gamma=2.4)	図2. 図1 を Gamut Mapping した結果(BT.2020色域, Gamma=2.4)



図3. 図1 を Gamut Mapping した結果(BT.709色域, Gamma=2.4)	図4. 図1 を 3x3 の Matrix で BT.2020 to BT.709 変換した結果(BT.709色域, Gamma=2.4)



図5. 図1と図2をアニメーションで切り替え(BT.2020色域, Gamma=2.4) 表示されない場合はココをクリック	図6. 図3と図4をアニメーションで切り替え(BT.709色域, Gamma=2.4) 　表示されない場合はココをクリック

図5 より BT.709 の Gamut Boundary の外側のデータのみが Mapping されていることが分かる。BT.709 の Gamut Boundary の内側のデータは変化しておらず意図した結果となっている。

加えて図6 より筆者にとっての従来手法である 3x3 の Matrix 変換と比較して、本手法では Hue のズレや色潰れが無くなっていることが分かる。

作成した 3DLUT は以下に置いた。興味があれば自由に使って頂いて問題ない（自分で言うのも何ですが DaVinci Resove の OpenFX の Color Space Transform を使った方が安心・安全だと考えます）。

drive.google.com

4. 理論

4.1. 概要

BT.2407 Annex2 の Gamut Mapping について簡単に説明する。

まず大前提として Gamut Mapping は CIELAB色空間で行う。IPT や CIECAM02 などを使わない理由は、BT.709 の色域外の「彩度が高い」領域で人間の視覚特性とマッチしない箇所があるためである。

さて、Gamut Mapping は以下の3つの要素から成り立っている。このうち Hue Mapping は冒頭で述べた理由により実装しない。

Lightness Mapping
Chroma Mapping
Hue Mapping <-- 今回は実装しない

ここから Lightness Mapping と Chroma Mapping の概要を説明する。この2つの Mapping は CIELAB色空間から算出できる Chroma-Lightness平面上で行う。この平面上で処理を行うことにより、Mapping の前後で CIELAB色空間での Hue は維持される。

BT.2407 Annex2 の Mapping 方法を説明する前に、そもそも Mapping の方法にはどういったものが考えられるか例を挙げておく。

図7. Chroma-Lightness 平面上での Gamut Mapping の例

図7 では BT.2020 色域のある点 (A) を BT.709 色域に Mapping している。図から読み取れるように Mapping の方法にはバリエーションが存在する。例えば (B) は Lightness の保持を最優先して、Chroma が大幅に減少したとしても Lightness を保持するようにしている。一方で (C) は Chroma を最優先して、Lightness が大幅に減少したとしても、Chroma を保持するようにしている。他にも Chroma-Lightness平面上で Mapping 前後のユークリッド距離を最小化する、などの方法も考えることができる。

さて、今回実装した BT.2407 Annex2 は (D) に示した方法となっている。(B), (C) の中間のような感じである。この座標がどのように算出されるのか以降で説明していく。

4.2. 詳細

先程の図7 の (D) に示した Mapping を実現するために、BT.2407 Annex2 では L_focal, C_focal と呼ばれる点を設定し、この点を基準として Gamut Mapping を行う。順を追って説明する。

4.2.1. L_focal の生成

L_focal を求めるには以下のステップを経由する。

BT.709 cusp と BT.2020 cusp を計算
上記の2つの cusp から L_cusp を計算
L_cusp の範囲を制限することで L_focal を生成

BT.709 cusp と BT.2020 cusp はそれぞれ Chroma-Lightness平面で最も Chroma が大きい点を意味し、L_cusp は BT.709 cusp と BT.2020 cusp を通る直線と Lightness軸との交点を意味する。文章での説明よりも動画の方がわかりやすいと思うので以下の図8 を参照して頂きたい。

図8. BT.709 cusp と BT.2020 cusp と L_cusp の関係

L_cusp が求まったら L_cusp の値を制限することで L_focal が生成できる（※4）。制限した結果を図9 に示す(図は筆者が書いたものではなく BT.2407 の Figure2-4 をコピペしたもの)。

※4 制限する明確な理由はなんとなくしか理解できてないので本記事では言及しない。ちなみに自分の実装では制限しないと破綻が起こった。

4.2.2. C_focal の算出

次に C_focal の算出方法について説明する。C_focal は BT.709 cusp と BT.2020 cusp を通る直線と Chroma軸との交点の絶対値である。したがって、以下のように2パターンが存在する。

パターン	図 (a) は BT.2407 の Figure A2-3 をコピペしたもの。(b) は Figure A2-3 を加工したもの
図10. パターン (a)
図11. パターン (b)

(a) は Chroma軸との交点が正の値となるケース、(b) は Chroma軸との交点が負の値となるケースである。いずれにせよ絶対値を取るため最終的な C_focal は正の数のみとなる。

4.2.3. L_focal, C_focal を基準とした Mapping 処理

L_focal, C_focal が決まった後は、この点を基準に Mapping を行う。はじめに具体例を図12 に示す。

図12 では BT.2020 色域の src point を BT.709 色域の dst point に Mapping している。図から分かるように、L_focal, C_focal を結ぶ直線よりも上側のデータは L_focal へ収束する直線上で、下側のデータは C_focal から発散する直線上で Mapping を行う。これが BT.2407 Annex2 の方式である（※5）。

Mapping先は BT.709 の Gamut Boundary である。この方法だと focal を基準とした直線上のデータは同一の Chroma, Lightness値に Mapping されるため色潰れが生じることを懸念するかもしれない。しかし色潰れが生じる可能性は極めて低い。実画像では Chroma が変化するとともに Lightness も変化するが、その変化が focal を基準とした直線に乗り続ける可能性が低いからである。この直線から外れると Mapping 後の値にも差が生じるため色潰れが生じることはない。

※5 この Mapping 方法が BT.2020 to BT.709 で適切な理由は分かっていない。個人的には L_focal を BT.709 cusp と BT.2020 cusp の中間の Lightness値にしても良いのでは？と思っている。

5. 実装

ここまでに述べてきた理論を実際に実装する。

5.1. 処理の大まかな流れ

初めに処理の大まかな流れを説明しておく。以下の図を参照して欲しい。

入力のRGB(Gamma=2.4)を Lab, LCH(Lightness, Chroma, Hue) に変換
該当する Hue の Chroma-Lightness平面をプロット
BT.709 cusp, BT.2020 cusp を算出
L_cusp, L_focal, C_focal を算出
入力の LCH から L_focal を使うのか C_focal を使うのか判別
判別した focal 基準で BT.709 の Gamut Boundary に Mapping
Mapping 後の LCH から RGB(Gamma=2.4)を算出する

項目	L_focal 基準の変換例	C_focal 基準の変換例
1
2
3, 4, 5
6
7

さて、上記の処理はそれなりにボリュームのある処理である。画像の全RGB値に対してバカ正直に処理すると重くなってしまう。

そこで今回の実装では Mapping処理を「入力のRGB値に依存しない処理」と「入力のRGB値に依存する処理」に分け、前者に関しては事前に計算しておき LUT化することにした。

具体的には以下のLUTを作成した。

Gamut Boundary 2DLUT
- BT.709, BT.2020 の Gamut Boundary 情報の入った 2DLUT
L_focal 1DlUT
- L_focal 情報の入った 1DLUT
C_focal 1DLUT
- C_focal 情報の入った 1DLUT
Chroma Mapping 2DLUT for L_focal
- 任意の Chroma-Lightness 平面での L_local からの角度 D_l のデータに対する L_focal からの距離 R_l の入った 2DLUT
Chroma Mapping 2DLUT for C_focal
- 任意の Chroma-Lightness 平面での C_local からの角度 D_c のデータに対する C_focal からの距離 R_c の入った 2DLUT

以降では、LUTの作成および使用方法を交えながら実装内容について説明していく。

5.2. BT.709, BT.2020 の Gamut Boundary を算出

基本的な理論は前回のブログで記した通り。ただし、この記事で書いた手法は遅すぎたため現在は別の手法を使用している。

trev16.hatenablog.com

5.3. BT.709 cusp, BT.2020 cusp の算出

4.2.1. で説明した条件を満たす点を算出するだけであり特記事項はない。5.2. で作った LUT から簡単に求めることが出来る。

5.4. L_focal LUT, C_focal LUT の算出

基本的には 5.3. で求めた BT.709 cusp, BT.2020 cusp を使って、4.2.1. , 4.2.2. で述べたとおりに機械的に求めるだけで良い。

ただし今回の実装では 5.2. で作成した Gamut Boundary の LUT の量子化誤差の影響で Hue に対して高周波成分が発生してしまった。真の値はガタついてないと推測されるため、LPF を適用して高周波成分を除去した。

また、C_focal に関しては LPF に加えて以下の2点の追加処理を実行した。

無限大になる値(※6)は付近の値から適当に補間
最大値を5000以下に制限(※7)

※6 これもLUT値の量子化誤差の影響。図では Zero Division Error として値を0にしてある
※7 C_focal の変化が一定値を超えた場合に後の線形補間計算での致命的な誤差が生じたため

結果を以下に示す。


図13. L_focal LUT (オレンジの線)	図14. C_focal LUT (オレンジの線)

計算した L_focal, C_focal は LUT化した。Hue を入力すると focal 値が出力される LUT とした。

5.5. focal を基準とした Mapping 先の算出

5.5.1. 方針

L_focal, C_focal が決まれば、次は Mapping 先の算出となる。繰り返しになるが図12 で示した通り Mapping 先は focal を基準とした直線と BT.709 の Gamut Boundary の交点である。したがって Mapping 処理の際は全ての入力RGB値（を変換した LCH値）に対して直線の数式を求めて BT.709 の Gamut Boundary との交点を計算する必要がある。しかし、これはかなり面倒である。

そこで今回は事前に focal から 1024本の直線を引き、 BT.709 の Gamut Boundary との交点を計算して LUT化しておくことにした。こうすることで、本番の計算時は単純な三角関数の計算とLUTの参照だけで Mapping 先の座標計算が可能となった。

以下で、生成した LUT の詳細を述べる。

5.5.2. 2D-LUT の作成

生成した LUT が概要は以下の通り。

LUT の入力は入力RGB値から計算した Hue と、focal からの角度 D とする
LUT の出力は focal からの距離 R とする
LUT の入力の D は最大値と最小値が存在し Hue によって値が変わる

上記から分かるように 2入力1出力の 2D-LUT である。また focal は L_focal, C_focal と2種類存在するので LUT も2種類作成した。各種パラメータを図にしたものを以下に示す。添字の "l" と "c" はそれぞれ L_focal用、C_focal 用のパラメータであることを意味する。

説明	図
図15. L_focal 用 LUT の説明
図16. C_focal 用 LUT の説明

focal からの角度 D は図にあるとおり、最小値と最大値を設けてある。当初は以下の通り固定値とすることを検討していた。

L_focal の D の範囲: -pi/2 ～pi/2
C_focal の D の範囲: 0 ～ pi

しかし、この範囲だと全く参照されない LUTデータの割合が多く補間計算の精度が低下した。そこで最小値と最大値が Hue に応じて動的に変わる仕様とした。 focal からの角度 D の範囲の詳細は以下の通り。

D_l_min: L_focal と C_focal を結ぶ直線に対する角度
D_l_max: pi/2 固定
D_c_min: C_focal と BT.2020 cusp を結ぶ直線に対する角度(※8)
D_c_max: pi 固定

※8 C_focal と L_focal ではなく、C_focal と BT.2020 cusp を結ぶ直線にした理由は C_focal の変化量が大きい Hue での計算誤差を防ぐため。本記事では詳細は言及しない。あくまでも筆者の実装上の都合である。

5.5.3. 2D-LUT の使用

作成した 2D-LUT は次のようにして使用した。

入力RGB値を LCH値に変換
Chroma-Lightness 平面上で D_l, D_c を計算
L_focal 用の 2D-LUT と C_focal 用の 2D-LUT に Hue, D_l, D_c を入力して focal からの距離 R_l, R_c を求める
入力LCH値が L_focal と C_focal を結ぶ直線の上側か下側かを判別。適切な方の結果を最終結果として採用
簡単な三角関数で Mapping 先の LCH値を計算
元の入力 LCH値が BT.709 Gamut Boundary の内側だった場合は計算結果を捨てる。

手順で書いた通り、最初は L_focal を使うか C_focal を使うか判別しない。両方の LUTから距離 R を求めている。これは筆者の実装上の都合であり別に最初に判別しても構わない。同様に入力RGB値が BT.709 Gamut Boundary の内側だった場合でも処理は一通り行う。これも筆者の実装上の都合である。

手順 3.と 4. の様子を可視化した動画を以下に示す。上段の途中計算の段階では LUT の想定範囲外の角度のデータは異常値となっているが、最終的には左下のように想定範囲内のデータのみが採用されるため問題ない。

図17. Gamut Mapping の途中経過の可視化

5.6. ソースコード

参考までにソースコードを載せておく。以下のリンクの bt2407_gamut_mapping_for_rgb_linear が該当の処理である。

github.com

6. 検証

65x65x65 の RGBパッチに対して本処理を適用し、異常が無いことを確認した。確認用の動画を図18に示す。

図18. 65x65x65 の RGBパッチに対する変換結果

動画で矢印が書かれいている点は実際に Gamut Mapping が行われたサンプルである。矢印が書かれていない点は Gamut Mapping が行われなかったサンプルである。BT.709 の Gamut Boundary の外側のデータは BT.709 の Gamut Boundary に Mapping されていることが分かる。また Gamut Boundary の内側のデータは Mapping されていないことが分かる。この結果から本実装が正しく行われていると判断する。

7. 考察というか感想

半年以上もかかってしまったがようやく Gamut Mapping を一つ実装できた。何も分かっていない素人の状態から少しは知識が身に付いたと思う。まだまだやりたいことが沢山あるので1つずつ着実にこなして行きたい。

今後の展望としては以下の点がある。

BT2446 と組み合わせた HDR10 to BT.709 変換の3DLUT の作成
- これは絶対にやる。すぐやる。たぶん1ヶ月以内にできる。
既に世の中に存在する既存のツールとの比較
- YouTube や DaVinci Resolve の Gamut Mapping との差異を確認したい
- 実は DaVinci Resolve との比較は軽く行っている。かなり近い変換結果になった。詳しく調べて結果をまとめたい。
任意の色域への拡張(DCI-P3 to BT.709 も実現できるようにする、など)
- 個人的には L_focal の値をもう少し制限すれば汎用化が可能だと推測している
- アイデアはあるので色々と試したい
高速化の検討
- 処理の後半で使用した 2D-LUT まわりの計算は LUTの生成方法も含めて簡易化をしたい

8. 参考文献

Report ITU-R BT.2407, "Colour gamut conversion from Recommendation ITU-R BT.2020 to Recommendation ITU-R BT.709", https://www.itu.int/pub/R-REP-BT.2407

toruのブログ

月1くらいで記事書きたい。

Report ITU-R BT.2407 の Annex2 を実装してみた