高动态范围成像

    自然景物的辐照范围很广，同样人眼也可以感知到很大的辐照范围下的图像，然而摄像机成像很容易饱和或者曝光不足，所以需要根据拍摄场景来设置合理的快门时间（或者光圈，ISO等参数）。

    尽管通过设置快门时间（或者光圈，ISO等参数）可以得到一幅比较理想的图像，但仍然可能存在曝光不足或曝光过度的现象。

    高动态范围成像（High Dynamic Range）通过

    1）拍摄一组不同曝光下的图像序列，根据已知曝光时间来估计相机的辐射响应函数，

    2）然后将所有像素点映射为辐照值，即得到高动态范围图像。

   但是由于显示设备都是以256阶强度显示图像，所以需要对高动态范围图像进行色调映射以适应显示设备，色调映射后的图像应该比不同曝光图像序列中任意图像所观察的动态范围大。

1 图像对齐

    在拍摄的一组不同曝光图像序列中，不同图像间不可避免存在一些位移，如果位移仅有相机抖动引起，可以使用全局平移参数进行描述。由于拍摄内容高度一致，所以可以使用一种简单方法估计平移，具体如下：

    1）求每帧图像的中值阈值图（Median Threshold Bitmap），由于曝光改变并不改变像素值大小关系，所以每帧图像的中值阈值图应该具有相同的模式，但可能存在一些位移；

    2）以图像序列中间帧作为参考，使用图像金字塔逐次求解每一帧求关于中间帧的全局平移参数；

    opencv 提供了类 AlignMTB 实现了这种情况下的图像对其，其计算全局平移主要代码如下：

// 创建图像金字塔
std::vector<Mat> pyr0;
std::vector<Mat> pyr1;
buildPyr(img0, pyr0, maxlevel);
buildPyr(img1, pyr1, maxlevel);

// 初始平移为0
Point shift(0, 0); 

// 自顶而下计算平移量
for(int level = maxlevel; level >= 0; level--) {
    shift *= 2;
    Mat tb1, tb2, eb1, eb2;
    computeBitmaps(pyr0[level], tb1, eb1);
    computeBitmaps(pyr1[level], tb2, eb2);

    int min_err = (int)pyr0[level].total();
    Point new_shift(shift);
    for(int i = -1; i <= 1; i++) {
        for(int j = -1; j <= 1; j++) {
             Point test_shift = shift + Point(i, j);
             Mat shifted_tb2, shifted_eb2, diff;
             shiftMat(tb2, shifted_tb2, test_shift);
             shiftMat(eb2, shifted_eb2, test_shift);
             bitwise_xor(tb1, shifted_tb2, diff);
             bitwise_and(diff, eb1, diff);
             bitwise_and(diff, shifted_eb2, diff);
             int err = countNonZero(diff);  // 非零元越少表示重合程度越高
             if(err < min_err) {
                 new_shift = test_shift;
                 min_err = err;
               }
            }
        }
     shift = new_shift;
  }

    如果图像位移由于拍摄场景内部的某些运动引起，则需要描述局部平移，这可以使用稠密光流等方式实现。在光流计算中，一般要求两帧间图像亮度变化较小，一种解决方案是在直方图均衡化后的图像序列上计算光流。

    结合全局平移与局部光流，可以对图像进行精确对齐。在HDR应用中，图像序列均在极短时间内获取，一般情况只需要计算全局平移即可。

2 求相机辐照响应函数

    当成像物体辐照强度不同时，经过相机后会产生不同的像素值，该系统可以用下图描述

    

    已知像素值，同时知道不同图像的快门时间，如何估计辐照响应函数呢？大致思路如下：

    对某一固定点 i，假设其辐照强度为 ，不同帧的快门时间为 ，不同帧上固定点 i 的像素值为 ，可建立关系 ，

    求反函数并取 2 的对数得 ，令 ， 即为需要求解的辐照响应函数（反函数）。

     即建立起一个最优求解方案，其中， 为未知函数， 为未知辐照强度。

    以上为相机辐照强度函数的基本求解思路，在实际实现中，能量函数上添加了一个平滑项以确保曲线平滑性。

    对于像素点的选取，总是利用亮度变化平缓区域，这可以避免梯度变化附近的光晕影响。

    opencv 提供了类 CalibrateDebevec 实现了相机辐照强度函数求解，如果图像为彩色图像，则得到3个不同通道的相机辐照函数。

    createCalibrateDebevec(int samples = 70, float lambda = 10.0f, bool random = false) 函数创建一个 CalibrateDebevec 对象，

    samples 使用足够采样点以获得一个超定的线性方程，lambda 为平滑像控制系数，random 决定采样点来自随机选取还是一个固定的矩形区域。

3 合成辐照图（radiance map）

    理想情况下，根据相机辐照响应函数，使用任意帧的非饱和点就可以求出任意点的辐照值。但真实图像中存在噪声干扰，可以使用加权方式求得该点的辐照值。

    

    观察以上辐照响应函数，在较暗区域或者饱和区域时，其梯度值趋近无穷大（注意，这里自变量为 Z），

    所以可以使用权重  进行加权处理， 。

    

    在每个通道上应用不同的相机辐照响应函数即可求得每个点上的辐照值，即构成了 HDR 图像。

    可使用 opencv 提供的 MergeDebevec 类合成辐照图，HDR 图像可以使用不同的格式存储，包括 .ppm，.hdr，.exr，每种格式对辐照值的表达上限不一致。

4 色调映射

    当获得 HDR 图像后，我们无法在显示器上直接显示结果，所以需要将 HDR 图像映射为 8 位深度图像。具体方法如下：

    1）将图像拆分为亮度通道与色彩通道；

    2）对亮度通道取对数 ；

    3）在亮度对数函数上应用低通滤波 ；

    4）亮度对数图像高频成分可表示为 ；

    5）压缩低频成分 ；

    6）形成新亮度对数函数 ；

    7）指数还原 ，即得到色调映射后结果；

    以上过程中低通滤波器如使用线性滤波（如高斯模糊）可能在图像边缘处产生光晕，这是由于线性滤波器不保留边界，可以使用双边滤波器改善光晕。

    双边滤波器不仅考虑像素间的位置信息，同时考虑像素间的亮度信息，这就较好的保留了边缘信息，从而避免边缘附近的光晕现象。

     

    opencv 在 photo 模块中提供了 HDR 相关实现，Tonemap 仅使用 gamma 调节，TonemapDrago 与 TonemapReinhard 引入了亮度压缩，同时可对色彩进行调整。

    以下给出示例图像：

        

       

参考资料   Computer Vision: Algorithms and Applications   Richard Szeliski

                  Recovering High Dynamic Range Radiance Maps from Photographs   Paul E. Debevec &Jitendra Malik