この記事は，P&D - Planning and Development - Advent Calendar 2017の21日目の記事です．
最近あまりP&Dに対してコミットできておらず幽霊部員となっていましたが，@3000manJPY 君からお誘い頂いたので記事を書いてみます．
投稿に不慣れなもので途中で画像のアップロードの制限がかかってしまい，本文内に作業途中の画像がいくつかあります．
後日差し替えますのでご了承ください．
2018/04/17 画像差し替え済
紹介する内容は，コンピュータビジョンにおける有名な技術の1つであるHigh-Dynamic-Range(HDR) Imagingについてです．
最近では当たり前のようにスマホに搭載されるようになったHDR撮影ですが，このHDRがどういう技術なのかを，実践を踏まえて説明することを目的とします．
この記事を書くに当たって初めて理解したところもあり，内容が正確でない箇所もあるかと思います．
もしも見つけた方がいらっしゃいましたら，ご指摘いただけると幸いです．
腑に落ちない点については，積極的に議論していけたらと思っています(優しくお願いします）．
実装環境
- Windwos 10
- Visualstudio 2015
- OpenCV 3.2.0
はじめに
私たちの住んでいる世界では，星空から直射日光までおよそ10桁以上の輝度レンジがあります．
その輝度レンジを人の目は限られた範囲でしか認識できませんが，瞳孔を閉じたり開いたりすることで，光を調整しより明るいシーンや暗いシーン対して適応しています（明暗順応）．
イメージは以下の図の通りです．
これと同様にカメラも，輝度レンジを認識する機械と考えることができます．
シャッタースピードやしぼりを調整することで露出を決定し，認識するシーンの輝度レンジを決めることも可能です（詳しくは次節で）．
ただし，人間よりは優れていないため狭い範囲(Low-Dynamic-Range(LDR))でしか認識するができません．
これらから分かるように，人間はHDRを識別できますが，カメラはLDRしか識別できません．
HDR輝度マップの作成
一方で，下図のように露出を変えた撮影画像をうまく使えば，それぞれの画像が得た輝度レンジを組み合わせることで，広範囲のレンジを獲得できるのではないか？と考えられると思います．
これが，HDR Imagingの根幹の考え方で，まずはこの輝度マップを作成すること目指します．
アイデアとしては単純で，カメラを固定して，シャッタースピードを変えた画像を複数枚撮影すれば，各画像の輝度情報を基に輝度マップを作成することができます．
（カメラを固定しないと各ピクセルの輝度マップが正しく算出できないので，注意が必要です．）
今回は２つのシーンに対して，輝度マップを作成しました．
１つ目は，ある部屋の窓辺で被写体を撮影したシーン，２つ目は夜景を撮影したシーンです．
- 夜景を撮影(Wikipediaから取得)
HDR Imagingの実装
各ピクセル(x,y)における，輝度マップを算出しています．
このとき，白とびや黒つぶれは輝度マップの作成に悪影響を与えているので除いた上で，作成しています．
注意する点としては，画素値÷露光時間を特定の条件で一定にしておく必要があることです．
実装
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ビット画像にすることができます．
トーンマッピングにはいくつか方法がありますが，今回は，単純な方法とOpenCVに実装されているものを試します．
実装と結果
単純な手法
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);
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);
OpenCVに実装されている手法
cv::Ptr<cv::TonemapReinhard> tonemap = cv::createTonemapReinhard(1.5, 0, 0, 0); // 適当なパラメータ
tonemap->process(hdr_image, ldr_image);
hdr_image *= 255;
実用上の問題点
ここまで話した内容は，画素値とシャッタースピードは比例の関係にある前提で扱ってきましたが，実際は違います．
これは一般的にカメラの撮像素子に電荷と出力される画素値が比例しないためです．
したがって各画素値について、非線形マッピング（応答関数と呼ばれる）が行われます．
実際にどのよう画素値と露光時間が対応しているかを知っていれば、シーンのダイナミックレンジを正しく再現することができます。
しかし，実際問題そんなこと滅多にありません．
そこで，この非線形の応答関数を推定し，線形に戻してから輝度マップを作成する手法を紹介します．
応答関数の推定と適応
今回使用したのは，Recovering High Dynamic Range Radiance Maps from Photographsで提案された手法です[Debevec,1997]．
まず，シャッタースピードがを変えた$j$番目の，ある画像の$i$番目の画素値$Z_ij$が輝度マップの$i$番目と$j$番目のシャッタースピードの掛け算を入力とした関数$f()$で与えられたとすると，
Z_ij = f(E_i\Delta t_j)
という式で定義できるとします．
このとき，未知なのは関数$f()$と$E_i$なのでこれらを求めます．
このままだと，直接求めることができないため，両辺で$log$をとって式変形をすると
g(Z_ij) = \ln{E_i} + \ln{\Delta t_j}
となります．
この考えもとに，既知の画素値とシャッタースピードを合わせて，以下の優決定性の線形方程式を解くことで，それぞれを求めます．
\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に実装されたものを使用しました．
実装
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を使用
左の画像では，配置されている物体が暗くならずに背景の雲も詳細に再現できている事が分かります．
また，右の画像でも同様に応答関数を考慮せずに行なった場合に比べると，アーチファクトが少なく良好な結果が得られています．
若干画像の色味が違っていますが，これはアルゴリズムによるもので，適切なパラメータを指定すれば改善されると思います（申し訳ありませんが，試せていません）．
考察
HDR Imagingは輝度マップを作成することが目的です．
輝度マップを作成することができれば，その応用の１つとしてトーンマッピングを用いた，画像表示が可能となります．
今回はパラメータ調整をしていないため，あまりうまく可視化ができたと言えないかもしれませんが，うまく調整することで適切なトーンマッピングを行う事ができ，人間の感覚として良いと思う画像になると思います．
もっと手軽に白とびや黒とびを軽減したいなら，Exposure Fusionを使ってみるのもありかもしれません．
これは，入力に露光時間を必要とせず明るさの変化するシーンの画像を与えれば，輝度マップを作成することなく画像を生成します．
cv::Ptr<cv::MergeMertens> merge_mertens = cv::createMergeMertens();
merge_mertens->process(images, fusion);
fusion*= 255;
Exposure Fusionの方が，簡単で露光時間の入力も必要ないので，とても合理的だと言えるでしょう．
しかし欠点もあります．
もし，好みのシーン画像を得られなかった場合には，微調整できるパラメータ等はありません．
その点，HDR Imagingは輝度マップを獲得しているので，トーンマッピングしだいてあらゆる可視化の方法を検討できます．
最後に
この記事では，HDR Imagingについて紹介しました．
ここで紹介したのはHDR Imaging技術の一部です．
この他にも輝度マップの作成や，トーンマッピングの手法はいくつかあります．
もし，また機会があればそれらも紹介できればと思っています．
参考
High Dynamic Range (HDR) Imaging using OpenCV
Recovering High Dynamic Range Radiance Maps from Photographs
High dynamic range imaging and tonemapping
おまけ
様々な露光時間で撮影した写真が少しでもずれると，上記の手法はうまく動作しません．
そこで，位置ズレを補正してくれる，Image Alignmentという技術を利用します．
実装
// 位置調整
cv::Ptr<cv::AlignMTB> alignMTB = cv::createAlignMTB();
alignMTB->process(images, images);
結果
