[CV]High-Dynamic-Range(HDR) Imagingについて

この記事は，[P&D - Planning and Development - Advent Calendar 2017](https://qiita.com/advent-calendar/2017/planningdev)の21日目の記事です．

最近あまり[P&D](https://www.planningdev.com/)に対してコミットできておらず幽霊部員となっていましたが，@3000manJPY 君からお誘い頂いたので記事を書いてみます．

~~投稿に不慣れなもので途中で画像のアップロードの制限がかかってしまい，本文内に作業途中の画像がいくつかあります．
後日差し替えますのでご了承ください．~~

2018/04/17 画像差し替え済


---

紹介する内容は，コンピュータビジョンにおける有名な技術の1つである**High-Dynamic-Range(HDR) Imaging**についてです．
最近では当たり前のようにスマホに搭載されるようになった**HDR**撮影ですが，この**HDR**がどういう技術なのかを，実践を踏まえて説明することを目的とします．

この記事を書くに当たって初めて理解したところもあり，内容が正確でない箇所もあるかと思います．
もしも見つけた方がいらっしゃいましたら，ご指摘いただけると幸いです．
腑に落ちない点については，積極的に議論していけたらと思っています(優しくお願いします）．

---

**実装環境**

+ Windwos 10
+ Visualstudio 2015
+ OpenCV 3.2.0

#はじめに
私たちの住んでいる世界では，星空から直射日光までおよそ**10桁以上**の**輝度レンジ**があります．
その輝度レンジを人の目は限られた範囲でしか認識できませんが，瞳孔を閉じたり開いたりすることで，光を調整しより明るいシーンや暗いシーン対して適応しています（明暗順応）．
イメージは以下の図の通りです．

![qiita_0.png](https://qiita-image-store.s3.amazonaws.com/0/211270/4caaec91-db55-aa7c-5d0c-07d6850c4fc4.png)

これと同様に**カメラ**も，輝度レンジを認識する機械と考えることができます．
シャッタースピードやしぼりを調整することで露出を決定し，認識するシーンの輝度レンジを決めることも可能です（詳しくは次節で）．
ただし，人間よりは優れていないため狭い範囲(**Low-Dynamic-Range(LDR)**)でしか認識するができません．

![qiita_1.png](https://qiita-image-store.s3.amazonaws.com/0/211270/a9125c4b-42e0-7733-a28e-f8b63a020dd5.png)

これらから分かるように，人間は**HDR**を識別できますが，カメラは**LDR**しか識別できません．

#HDR輝度マップの作成
一方で，下図のように露出を変えた撮影画像をうまく使えば，それぞれの画像が得た輝度レンジを組み合わせることで，広範囲のレンジを獲得できるのではないか？と考えられると思います．
これが，**HDR Imaging**の根幹の考え方で，まずはこの輝度マップを作成すること目指します．

アイデアとしては単純で，カメラを固定して，シャッタースピードを変えた画像を複数枚撮影すれば，各画像の輝度情報を基に**輝度マップ**を作成することができます．
（カメラを固定しないと各ピクセルの輝度マップが正しく算出できないので，注意が必要です．）




![qiita_2.png](https://qiita-image-store.s3.amazonaws.com/0/211270/03e8f845-9ea1-c864-f096-237571514a54.png)

今回は２つのシーンに対して，輝度マップを作成しました．
１つ目は，ある部屋の窓辺で被写体を撮影したシーン，２つ目は夜景を撮影したシーンです．

+ 夜景を撮影([Wikipediaから取得](https://en.wikipedia.org/wiki/High-dynamic-range_imaging))
<img alt="qiita_4.png" src="https://qiita-image-store.s3.amazonaws.com/0/211270/322ea40b-0b57-54eb-8ac7-41cf9be9b3a8.png">

##HDR Imagingの実装
各ピクセル(x,y)における，輝度マップを算出しています．
このとき，白とびや黒つぶれは輝度マップの作成に悪影響を与えているので除いた上で，作成しています．
注意する点としては，画素値÷露光時間を特定の条件で一定にしておく必要があることです．

##実装

```C++:HDR-Imaging
cv::Mat createHDR_naive(const vector<cv::Mat>& images, const cv::Mat& times) {
	cv::Mat image;

	auto xmin = 5, xmax = 250; // 白とび，黒つぶれを定義

	image = cv::Mat::zeros(images[0].size(), CV_32FC3);
	
	for (int y = 0; y < image.size().height; y++) {
		for (int x = 0; x < image.size().width; x++) {
			cv::Vec3i non_num = 0; cv::Vec3f non_sum = 0.0f;
			for (int i = 0; i < static_cast<int>(images.size()); i++) {	
				auto pixel = images[i].at<cv::Vec3b>(y, x);			
				for (int j = 0; j < pixel.rows; j++) {
					if (pixel[j] > xmin && pixel[j] < xmax ) {
						non_sum[j] += pixel[j] / times.at<float>(i); // 露光時間で割って一定に
						non_num[j]++;
					}
				}
			}
			for (int j = 0; j < non_num.rows; j++) {
				if (non_num[j] > 0)
					image.at<cv::Vec3f>(y, x)[j] = non_sum[j] / non_num[j];
			}
		}
	}

	return std::move(image);
}
```


(もっとスッキリ書きたかったのですが，~~グレースケール画像を想定したものを拡張したため，ネストの深い煩雑なコードになった事を~~お許しください．)

この関数によって，返されるのが輝度マップとなります．

#HDRの可視化
輝度マップが得られたからといって、ディスプレイ等の出力装置に直接表示することはできません．
なぜなら，今回作成した輝度マップは32ビットで表現されているためです．
そこで，**トーンマッピング**を行う事で，ディスプレイに表示可能な**8ビット画像**にすることができます．

![qiita_21.png](https://qiita-image-store.s3.amazonaws.com/0/211270/0908e992-6278-76f2-0838-5a02a1d26dd7.png)

トーンマッピングにはいくつか方法がありますが，今回は，単純な方法とOpenCVに実装されているものを試します．

##実装と結果
###単純な手法
```C++:平均値が128になるような変換	
auto mean = cv::mean( HDR_image);
hdr_image.forEach<cv::Vec3f>([mean](cv::Vec3f &x, const int *position) -> void {
	for (int i = 0; i < x.rows; i++)
		x[i] *= 128.0f / mean[i];
});
HDR_image.convertTo(HDR_image, CV_8UC3);
```
![mean.png](https://qiita-image-store.s3.amazonaws.com/0/211270/296eb527-3641-62f3-2b47-3452e7f08776.png)



```C++:ガンマ変換
double min, max;
auto gamma = 2.2;
cv::minMaxLoc(HDR_image, &min, &max);
HDR_image = (HDR_image - min) / (max - min);
HDR_image.forEach<cv::Vec3f>([](cv::Vec3f &p, const int *position) -> void {
	for(auto &x : p.val)
		x = cv::pow(x, gamma);
});
HDR_image.convertTo(HDR_image *= 255.0f, CV_8UC3); 
```
![gamma.png](https://qiita-image-store.s3.amazonaws.com/0/211270/9c2f8f5b-f344-7709-2d42-ff09298d7e44.png)

###OpenCVに実装されている手法

```C++:Reinhard
cv::Ptr<cv::TonemapReinhard> tonemap = cv::createTonemapReinhard(1.5, 0, 0, 0); // 適当なパラメータ
tonemap->process(hdr_image, ldr_image);
hdr_image *= 255;
```
![re.png](https://qiita-image-store.s3.amazonaws.com/0/211270/f4d61e62-efcc-5c62-08fd-e11923edcd4c.png)



#実用上の問題点
ここまで話した内容は，画素値とシャッタースピードは比例の関係にある前提で扱ってきましたが，実際は違います．
これは一般的にカメラの撮像素子に電荷と出力される画素値が比例しないためです．

<img alt="image" src="https://qiita-image-store.s3.amazonaws.com/0/211270/38886bab-813a-c1aa-367d-2434a4b7217a.png">

したがって各画素値について、非線形マッピング（**応答関数**と呼ばれる）が行われます．
実際にどのよう画素値と露光時間が対応しているかを知っていれば、シーンのダイナミックレンジを正しく再現することができます。
しかし，実際問題そんなこと滅多にありません．

そこで，この非線形の**応答関数**を推定し，線形に戻してから**輝度マップ**を作成する手法を紹介します．

#応答関数の推定と適応
今回使用したのは，[Recovering High Dynamic Range Radiance Maps from Photographs](http://www.pauldebevec.com/Research/HDR/debevec-siggraph97.pdf)で提案された手法です[Debevec,1997]．
まず，シャッタースピードがを変えた$j$番目の，ある画像の$i$番目の画素値$Z_ij$が輝度マップの$i$番目と$j$番目のシャッタースピードの掛け算を入力とした関数$f()$で与えられたとすると，

```math
Z_ij = f(E_i\Delta t_j)
```

という式で定義できるとします．
このとき，未知なのは関数$f()$と$E_i$なのでこれらを求めます．
このままだと，直接求めることができないため，両辺で$log$をとって式変形をすると

```math
g(Z_ij) = \ln{E_i} + \ln{\Delta t_j}
```

となります．
この考えもとに，既知の画素値とシャッタースピードを合わせて，以下の優決定性の線形方程式を解くことで，それぞれを求めます．

```math
\sum_{i=1}^N \sum_{j=1}^P [\ln{E_i} + \ln{\Delta t_j}-g(Z_{ij})]^2+\lambda \sum_{z=Z_{min}}^{Z_{max}} g^{''}(z)^2
```

コードは論文の方に記載されているので，興味のある方はそちらを参考にしてください．

自前で実装してみたのですが，パフォーマンスがかなり悪かったので，OpenCVに実装されたものを使用しました．

###実装
```C++

vector<cv::Mat> images; // 画像群
cv::Mat times; // シャッタースピード

~ 中略 ~

// 応答関数の推定
cv::Mat response;
cv::Ptr<cv::CalibrateDebevec> calibrate = cv::createCalibrateDebevec();
calibrateDebevec->process(images, response, times); 

// 応答関数を考慮してHDR画像を作成
cv::Mat hdr_image;
cv::Ptr<cv::MergeDebevec> mergeDebevec = cv::createMergeDebevec();
mergeDebevec->process(images, hdr_image, times, response);
```

最終的には`hdr_image`に，HDR輝度マップがが格納されています．

###結果
**Reinhardを使用**

![debe.png](https://qiita-image-store.s3.amazonaws.com/0/211270/8e0e001d-ea5d-933a-b6c0-a6c11bbafd10.png)

左の画像では，配置されている物体が暗くならずに背景の雲も詳細に再現できている事が分かります．
また，右の画像でも同様に応答関数を考慮せずに行なった場合に比べると，アーチファクトが少なく良好な結果が得られています．
若干画像の色味が違っていますが，これはアルゴリズムによるもので，適切なパラメータを指定すれば改善されると思います（申し訳ありませんが，試せていません）．

#考察
**HDR Imaging**は輝度マップを作成することが目的です．
輝度マップを作成することができれば，その応用の１つとしてトーンマッピングを用いた，画像表示が可能となります．
今回はパラメータ調整をしていないため，あまりうまく可視化ができたと言えないかもしれませんが，うまく調整することで適切なトーンマッピングを行う事ができ，人間の感覚として良いと思う画像になると思います．

もっと手軽に白とびや黒とびを軽減したいなら，**Exposure Fusion**を使ってみるのもありかもしれません．
これは，入力に露光時間を必要とせず明るさの変化するシーンの画像を与えれば，輝度マップを作成することなく画像を生成します．

```C++:Expose-Fusion
cv::Ptr<cv::MergeMertens> merge_mertens = cv::createMergeMertens();
merge_mertens->process(images, fusion);
fusion*= 255;
```
![fusion.png](https://qiita-image-store.s3.amazonaws.com/0/211270/0d6df91d-8100-078e-e77c-ecacfb967bc8.png)


**Exposure Fusion**の方が，簡単で露光時間の入力も必要ないので，とても合理的だと言えるでしょう．
しかし欠点もあります．
もし，好みのシーン画像を得られなかった場合には，微調整できるパラメータ等はありません．
その点，**HDR Imaging**は輝度マップを獲得しているので，トーンマッピングしだいてあらゆる可視化の方法を検討できます．


#最後に
この記事では，**HDR Imaging**について紹介しました．
ここで紹介したのは**HDR Imaging**技術の一部です．
この他にも**輝度マップ**の作成や，**トーンマッピング**の手法はいくつかあります．

もし，また機会があればそれらも紹介できればと思っています．

#参考
[High Dynamic Range (HDR) Imaging using OpenCV](https://www.learnopencv.com/high-dynamic-range-hdr-imaging-using-opencv-cpp-python/)
[Recovering High Dynamic Range Radiance Maps from Photographs](http://www.pauldebevec.com/Research/HDR/)
[High dynamic range imaging and tonemapping](http://cybertron.cg.tu-berlin.de/eitz/hdr/index.html)

#おまけ
様々な露光時間で撮影した写真が少しでもずれると，上記の手法はうまく動作しません．
そこで，位置ズレを補正してくれる，**Image Alignment**という技術を利用します．

##実装
```C++
// 位置調整
cv::Ptr<cv::AlignMTB> alignMTB = cv::createAlignMTB();
alignMTB->process(images, images);
```
##結果
<img alt="" src="https://qiita-image-store.s3.amazonaws.com/0/211270/915b3377-68ad-3f75-f3cc-0891f10bb294.png">