ISP-不同的调光方法（ToneMapping）线性、直方图、限制自适应直方图

ISP-不同的调光方法（Tone Mapping）线性、直方图、限制自适应直方图

调光，更一般的称为Tone Mapping，用于将颜色从原始色调（通常是高动态范围，HDR）映射到目标色调（通常是低动态范围，LDR），映射的结果通过介质进行显示，在人眼的视觉特性的作用下，达到尽可能复原原始场景的效果。Tone mapping的效果取决于从原始色调到目标色调的映射关系，这一映射关系可以是简单的缩放或截断，也可以是对摄影中已有关系的逼近。
如果直接对数值进行截断，毫无疑问会丢失场景中明亮部分的细节，因为它们都被映射成为白色（或者其它处于范围上限的颜色）。简单的缩放，例如将场景中的所有颜色都除以场景中颜色的最大值。当场景中一小部分很亮，其它部分相对很暗时，简单的缩放会导致暗部的细节丢失。
简而言之调光\Tone Mapping在ISP的作用是将相机传感器捕捉到的14位的图像数据压缩到显示器能够呈现的8位
常见的调光方法有：
1、线性调光
2、直方图调光
3、限制自适应直方图调光
4、混合调光

1、线性调光

线性调光的思路很简单，就是将原始的像素区间线性的压缩至目标的像素区间，比如将14位0-16383的像素区间，压缩至0-256的像素区间，代码如下

def lineToMapping(y16, rate=0.1):
    y16_flatten = y16.flatten()
    y16_sorted = np.sort(y16_flatten)
    rate_num = int(len(y16_sorted) * rate)
    y16_max = y16_sorted[len(y16_sorted) - rate_num]
    y16_min = y16_sorted[rate_num]
    y16_rate = (y16 - y16_min) / (y16_max - y16_min)  #这里是源数据的动态范围的选择
    #y16_rate = y16  / 65535  // 如果确保源数据是0-65535并且线性映射0-255就采用这行代码
    y8 = np.clip(y16_rate * 255, 0, 255).astype(np.uint8)
    return y8

代码非常简单，rate就是抛点率，上下同时抛弃多少百分比的点。为什么要抛点呢？是因为防止14位数据中过亮或者过暗的点会占据较大的线性空间，调光后的图像过亮或者过暗对比度很差。线性调光虽然很简单，抛点也能解决部分对比度较低问题，但经线性调光后的图像的对比度往往较低，效果不是很理想。
如何使一幅图像的像素分布更加的均匀呢？对，直方图均衡。

2、直方图调光

直方图均衡常用语对比度不高的图像，目的是将一幅图像的直方图从一个分布尽量拉成均匀分布。那么在这个过程中，我们也可以将一个14位的图像空间，既映射到了8位空间，同时还使得像素的分布更加均衡增加了图像的对比度。代码如下

def HEToMapping(y16):
    y16 = y16.astype(int)
    y16_flatten = y16.flatten()
    hist_ = np.zeros((1, 16383))
    # 求直方图统计
    for i in range(len(y16_flatten)):
        a = y16_flatten[i]
        hist_[0, a] += 1
    hist_1 = hist_ / hist_.sum()
    CDF_hist = np.cumsum(hist_1)
    hist_new = CDF_hist * 255
    hist_new = hist_new.astype(np.uint8)
    y16 = hist_new[y16]
    return y16.astype(np.uint8)

上述代码就实现了利用直方图均衡化，将14bit的数据映射到8bit空间，同时与线性调光的纯线性映射不同，直方图调光会使得全局的对比度更大，细节更加的突出。
但是直方图调光也存在它的缺点，那就是局部的对比度还是不够好，那怎么样进一步提高局部对比度呢？有图像处理经验的同学应该可以想到，局部往往对应的就是分块，也就是接下来的方法。

3、限制自适应直方图调光

如果仅仅只是分块利用直方图调光再将调光后的分块图像拼接到一起，局部对比度是有提升但是块与块之间对比度不同就会有着明显的界线，图像会出现分块情况，为了防止这种情况，就需要限制每一块的对比度，同时每个像素都对应着自己的块，但是这样会使得计算量很大，所以一般的做法是将图像分成n×n的块，其余的像素采用插值的方法计算其映射后的值。代码如下

def CLAHEToMapping(y16, r=40, step_=8):
    y16 = y16.astype(int)
    width = y16.shape[1]
    height = y16.shape[0]
    width_block = width // step_  # 列
    height_block = height // step_  # 行   512行 640列

    total = width_block * height_block

    # 保留直方图
    hist_ = np.zeros((step_ * step_, 16383), dtype=np.float64)
    # 保留累积直方图
    CDF_hist = np.zeros((step_ * step_, 16383), dtype=np.float64)
    for i in range(step_):
        for j in range(step_):
            y16_block = y16[i * height_block:(i + 1) * height_block, j * width_block:(j + 1) * width_block]
            y16_block = y16_block.flatten()
            num = i * step_ + j
            for ii in range(width_block * height_block):
                hist_[num, y16_block[ii]] += 1
            LIMIT = r * width_block * height_block // 16383
            sum = 0
            for ii in range(16383):
                if hist_[num, ii] >= LIMIT:
                    sum += hist_[num, ii] - LIMIT
                    hist_[num, ii] = LIMIT
            bonus = sum // 16383
            for ii in range(16383):
                hist_[num, ii] += bonus
            # 获得累积直方图
            for ii in range(16383):
                if ii == 0:
                    CDF_hist[num, ii] = hist_[num, ii] / total
                else:
                    CDF_hist[num, ii] = CDF_hist[num, ii - 1] + hist_[num, ii] / total

    for i in range(height):
        for j in range(width):
            # 四个角区域
            # 左上
            if (i <= height_block / 2) and (j <= width_block / 2):
                y16[i, j] = CDF_hist[0, y16[i, j]] * 255
            # 右上
            elif (i <= height_block / 2) and (j >= ((step_ - 1) * width_block + width_block / 2)):
                y16[i, j] = CDF_hist[step_ - 1, y16[i, j]] * 255
            # 左下
            elif (i >= ((step_ - 1) * height_block + height_block / 2)) and (j <= width_block / 2):
                y16[i, j] = CDF_hist[(step_ - 1) * step_, y16[i, j]] * 255
            # 右下
            elif (i >= ((step_ - 1) * height_block + height_block / 2)) and (
                    j >= ((step_ - 1) * width_block + width_block / 2)):
                y16[i, j] = CDF_hist[step_ * step_ - 1, y16[i, j]] * 255
            # 四边区域 线性插值
            # 左
            elif (i >= height_block / 2) and (i <= (step_ - 1) * height_block + height_block / 2) and (
                    j <= width_block / 2):
                num_j = 0
                num_i = (i - height_block / 2) // height_block
                num1 = num_i * step_ + num_j
                num2 = num1 + step_
                num1 = int(num1)
                num2 = int(num2)
                p = (i - (num_i * height_block + height_block / 2)) / (1.0 * height_block)
                q = 1 - p
                y16[i, j] = CDF_hist[num1, y16[i, j]] * 255 * q + CDF_hist[num2, y16[i, j]] * 255 * p
            # 右
            elif (i >= height_block / 2) and (i <= (step_ - 1) * height_block + height_block / 2) and (
                    j >= (step_ - 1) * width_block + width_block / 2):
                num_j = step_ - 1
                num_i = (i - height_block / 2) // height_block
                num1 = num_i * step_ + num_j
                num2 = num1 + step_
                num1 = int(num1)
                num2 = int(num2)
                p = (i - (num_i * height_block + height_block / 2)) / (1.0 * height_block)
                q = 1 - p
                y16[i, j] = CDF_hist[num1, y16[i, j]] * 255 * q + CDF_hist[num2, y16[i, j]] * 255 * p
            # 上
            elif (i <= height_block / 2) and (j >= width_block / 2) and (
                    j <= (step_ - 1) * width_block + width_block / 2):
                num_j = (j - width_block / 2) // width_block
                num_i = 0
                num1 = num_j
                num2 = num1 + 1
                num1 = int(num1)
                num2 = int(num2)
                p = (j - (num_j * width_block+width_block/2)) / (1.0 * width_block)
                q = 1 - p
                y16[i, j] = CDF_hist[num1, y16[i, j]] * 255 * q + CDF_hist[num2, y16[i, j]] * 255 * p
            # 下
            elif (i >= (step_ - 1) * height_block + height_block / 2) and (j >= width_block / 2) and (
                    j <= (step_ - 1) * width_block + width_block / 2):
                num_i = step_ - 1
                num_j = (j - width_block / 2) // width_block
                num1 = num_i * step_ + num_j
                num2 = num1 + 1
                num1 = int(num1)
                num2 = int(num2)
                p = (j - (num_j * width_block+width_block/2)) / (1.0 * width_block)
                q = 1 - p
                y16[i, j] = CDF_hist[num1, y16[i, j]] * 255 * q + CDF_hist[num2, y16[i, j]] * 255 * p
            else:
                # # 其余点 需要双线性插值
                num_i = (i - height_block / 2) // height_block
                num_j = (j - width_block / 2) // width_block
                num_i = int(num_i)
                num_j = int(num_j)
                num1 = num_i * step_ + num_j
                num2 = num1 + 1
                num3 = num1 + step_
                num4 = num2 + step_
                u = (j - (num_j * width_block + width_block / 2)) / (1.0 * width_block)
                v = (i - (num_i * height_block + height_block / 2)) / (1.0 * height_block)
                y16[i, j] = ((1 - u) * (1 - v) * CDF_hist[num1, y16[i, j]] * 255 +  u * (1-v) * CDF_hist[
                    num2, y16[i, j]] * 255 + v * (1-u) * CDF_hist[num3, y16[i, j]] * 255 + u * v * CDF_hist[
                                 num4, y16[i, j]] * 255)
                # # if y16[i, j] > 255:
                # #     y16[i, j] = 255
                # # elif y16[i, j] < 0:
                # #     y16[i, j] = 0
            y16[i, j] = y16[i, j] + (y16[i, j] << 8) + (y16[i, j] << 16)
    return y16.astype(np.uint8)

双线性插值原理见这里
限制自适应直方图调光可以很好的解决图像局部对比度的问题，但是可能会带来图像总体亮度偏暗，这时将线性调光的结果加权融合就可以很好的解决这个问题。

4、混合调光

如上所说就是将两种调光的y8按照加权的方式融合即可。

5、总结

1、上述的每种调光都有自己的优点，可以综合考虑。
2、还有一些更新的tone mapping的方法，可以看看相关的论文。