ISP流程--LSC详解

目录

1.1. 基本概念

 

梳理如下：

1.1.1. 影响因素：linearization

1.1.3. 影响因素：OB

1.1.4. 影响因素：AWB

1.1.5. 影响因素：lens和sensor

1.2. 算法设计

1.2.1. ISP pipeline

1.2.2. 矫正方法

 1.2.3. LSC代码

1.2.4. 联动方式

1.3. 调试方法

1.4. 测试方法

1.4.1. 客观测试

1.4.1.1. 测试设备

1.4.1.2. 测试标准

1.4.2. 主观测试

---

 

1.1. 基本概念

我们首先了解一下lens shading，lens shading是由于Lens（镜头）的物理结构原因引起的。由于Sensor边缘接受到的光线比sensor中心弱，从而造成了中心亮四周暗的现象（有的文档也称这种现象为“渐晕”）。

 

Luma Shading就是我们常说的Vignetting，形成原因可以分为大致两个原因：

其一，由镜头本身引入。镜头本身就是一个凸透镜，由于凸透镜原理，中心的感光必然比周边多；

其二，与光的特性有关。根据光的衰减特性，光的衰减与入射角呈cos关系，所以平行于光轴的光能量保留的最多，垂直于光轴入射的光线为0。这两个原因导致画面中间亮，边缘暗。

 

 

Color Shading的成因相对Luma Shading会更复杂一些。形成原因大致可以分为三部分：

其一：与光的特性有关。我们都知道光的色散现象，即白光通过三棱镜后会被分解为七色光。而产生这种现象的原因就是三棱镜对不同波长光线的折射不同，从而导致不同波长光线走过的光程不同。同样的由于镜头对不同光谱光线的折射程度不同，导致入射光线中不同波长的光线落在Sensor的不同位置，从而引起Color Shading。

其二：由红外滤光片引入。Sensor可以感知人眼感知不到的红外光和紫外光，因此需要使用另外的滤波片进行滤除，否则会导致红绿蓝像素点的亮度值与人眼观察到的亮度值存在较大的差异。红外截止滤波片位于镜头和图像传感器之间，为了消除这些投射到Sensor上不必要的光线，防止Sensor产生伪色/波纹，从而提高色彩还原性。红外滤波片主要分为干涉型、吸收型。红外滤光片为干涉型红外截止滤波片，即反射红外光。在可见光区域有较高的透过率，存在较低反射率，而在红外区域正好相反，反射率较高，透过率很低。但成角度拍摄照片时，红外光在IR（Infrared简称IR）膜上会有较大反射，经过多次反射后，被Sensor接收从而改变图像R通道的值，引起Color Shading。蓝玻璃则是吸收型红外截止滤波片，对红外光有很强的吸收作用，不存在很大的反射，能在一定程度上减轻shading和色差问题。

其三：与镜头相关。由Sensor上微透镜的CRA（Chief ray angle - 主光线角）与镜头的CRA不匹配导致。镜头的主光线角与传感器不匹配，会使传感器的像素出现在光检测区域周围，致使像素曝光不足，亮度不够。

 

备注：CRA是Chief Ray Angle的缩写，意思是主光角。从镜头的传感器一侧，可以聚焦到像素上的光线的最大角度被定义为一个参数，称为主光角(CRA)。对于主光角的一般性定义是：此角度处的像素响应降低为零度角像素响应(此时，此像素是垂直于光线的)的80%。

 

我们在挑选Lens的时候会有一个CRA的参数，在选择sensor的时候同样有一个CRA的参数，一般我们要求Lens的CRA曲线与senosr的CRA曲线完全匹配，如果实在不能匹配，我们也要求在同样的像高位置，CRA相差不能超过3度，而且最好是Lens的CRA比sensor的CRA小。否则将会导致成像照度或色彩问题。

 

梳理如下：

镜头阴影现象及其成因：

日常生活中，我们接触到的摄像头在成像过程中往往会受到很多因素的干扰，例如坏点、畸变、噪声等等，ISP的一个重要任务便是要消除这些因素对最终成像的干扰，获得高质量的输出图像。

 

镜头阴影也是这些常见干扰因素之一，顾名思义，它会在摄像头最终成像上造成一片阴影。ISP中通常会设置镜头阴影矫正模块来处理这个问题。

 

那么镜头阴影的具体形态以及成因到底是什么？下面让我们来介绍一下。

阴影的种类与成因：

镜头阴影（Lens Shading）通常指的是图像中亮度或颜色的分布不均。最常见的镜头阴影表现为图像的四个角比中心区域要暗，或者四周的颜色与中心不同。这种现象在光学摄影中非常普遍，通常分为两大类：亮度阴影（Luma Shading）和色彩阴影（Color Shading）。

1. Luma Shading

Luma Shading的形态如下图所示。它体现为图像中心区域较亮，而四周较暗。

 

这种阴影的形成主要与镜头的结构和光学特性有关，具体可以归因为两点。

第一，镜头本身机械结构导致，该原因造成的阴影称为机械阴影（Mechanical Shading）。由于图像边缘的点通常对应较大的入射光角度，这些大角度光线在镜头中传播的时候，会更容易被镜筒、镜片边缘拦截，这会使得最后成像的光束大大减小，因而导致图像边缘亮度相比起中心有所下降。如下图所示，原本光束的直径为D2，经过拦截后，最终抵达sensor的光束直径只剩D3。

 

第二，镜头本身光学特性导致，该原因造成的阴影称为自然阴影（Natural Shading）。即便镜头不存在机械阴影，最后成像的四周仍然可能比中心暗，这主要是因为镜头中心与边缘的聚光能力不同。具体来说，对于整个镜头，我们可将其视为一个凸透镜，由于凸透镜中心的聚光能力大于边缘，因而会导致sensor中心接收的光线强度大于四周。这种从中心向四周衰减的速率符合余弦四次方定律，可以通过调整透镜系统来改良这种现象。余弦四次方定律的公式如下：

其中，为中心光强，θ为斜射到成像面上的光线与水平轴的夹角。

2. Color Shading

Color Shading的形态如下图所示。它的特征为图像中心与四周颜色不一致，出现了偏色的问题。

 

这种偏色现象通常是由镜头内部的光学结构、传感器的光学特性等因素共同作用引起的。具体可以归因为三点。

第一，红外截止滤波片（IR-Cut Filter）的引入。在相机中，红外截止滤波片位于镜头与sensor之间，它的作用是去除不必要的红外光，以免最终RGB像素点的亮度值与人眼观察到的亮度值存在较大的差异，引入颜色偏差。常见的红外截止滤波片利用反射原理来消除红外光，然而，当我们成一定角度地拍摄照片时，红外光在滤光片上会有较大角度反射，经过多次反射后，最终可能又会被sensor接收到，进而引起图像偏色。

第二，光的色散现象。由于镜片对不同波长光线的折射率不同，会导致各种颜色的光shading程度轻重不同，进而导致色彩失真。

第三，sensor上microlens的CRA（主光线角）与镜头的CRA不匹配。CRA的差异会导致光线照射到传感器时产生位置偏差，可能会错误地折射到周边像素的位置上，导致像素之间出现干扰，进而影响图像的色彩和亮度。

1.1.1. 影响因素：linearization

一般来说画面局部偏色还可能与linearization相关，评估sensor时应关注sensor线性化的特性，部分DSP处理中也会在ISP pipeline的开始加入sensor linearization的矫正。

 

对于人眼来说，shading系数大于80%，人眼就不太看得出来暗角。补偿的过高反而会因为四角增益放的过大，导致四角噪声偏大，暗处四角偏红等（偏红原因与BLC模块有关）。目前手机影像ISP流程的LSC补偿一般都放到85%左右，这样可以做到主观上人眼没有感受到明显的暗角，四角的噪声也不至于过大，是一个权衡的选择。

 

目前手机平台主流LSC的做法是挑选出一批模组样本中最具有代表性的golden模组，分块统计出golden模组的亮度信息，然后把这些信息填入到效果代码里。每一颗模组的亮度信息也会烧录在自身的OTP里。这样针对每一颗模组来说，都可以通过比较自身的亮度信息和代码里面的golden信息，从而调节自身的补偿系数，达到模组效果一致的目的。

1.1.3. 影响因素：OB

部分国产sensor的性能较差，画面边缘的OB可能和中心不同（高温下尤其明显），如果OB仅整体扣除一个数值，反映到最终图像效果上可能会造成lens shading或者color shading。

1.1.4. 影响因素：AWB

AWB并不会造成color shading，但是AWB会放大color shading的影响。

1.1.5. 影响因素：lens和sensor

Shading的根本影响因素是lens和sensor，物料选型时应注意sensor、lens的CRA匹配，sensor靶面大小与lens匹配。

1.2. 算法设计

1.2.1. ISP pipeline

Shading一般在OB和DPC的后面。额外注意，如果3A的统计数据在shading之后获取，那么shanding会影响3A的统计数据。

1.2.2. 矫正方法

目前手机影像主流LSC的做法是拍摄一张亮度适中均匀的图，将图像分成mxn块，计算每一块的平均亮度，然后将这些亮度信息存储在OTP里面，平台ISP处理时会把这些亮度信息读出来，然后计算每一块需要补偿的增益，再通过cos4次方曲线或者其他方式拟合一个补偿的矩阵乘到原图像上。

 

Shading的矫正方法目前主流有两种：一种是同心圆法，一种网格法。

下面介绍一下网格法流程：

 

• 图像划分：将原始图像划分为多个小的网格区域。每个网格的大小由用户指定，这样可以根据具体需求调整网格的密度。
• 计算网格均值：对每个网格区域的像素值进行平均计算，得到每个网格的平均亮度值。
• 计算增益矩阵：选取某个基准网格（通常是中心网格）的平均亮度值，计算其他网格相对于基准网格的亮度增益值。
• 插值计算：对于网格之间的像素点，使用插值法计算其增益值，以保证校正后的图像平滑过渡。
• 应用增益：将计算得到的增益矩阵应用于原始图像，对每个像素的亮度进行调整，得到光照均衡校正后的图像。

通过以上步骤，网格矫正法能够有效地校正图像的光照不均匀性，提高图像的整体质量。接下来，我们将通过 Python 代码演示如何实现这一方法。

import cv2
import numpy as np
import matplotlib.pyplot as plt
 
 
def calculate_lsc_gain(image, side_num):
    height, width, _ = image.shape
    side_y = height // side_num
    side_x = width // side_num
 
    # Separate color channels
    image_r = image[:, :, 0]
    image_g = image[:, :, 1]
    image_b = image[:, :, 2]
 
    # Initialize matrices to store block means
    image_point_r = np.zeros((side_num + 1, side_num + 1))
    image_point_g = np.zeros((side_num + 1, side_num + 1))
    image_point_b = np.zeros((side_num + 1, side_num + 1))
 
    # Calculate block means
    for i in range(side_num + 1):
        for j in range(side_num + 1):
            x_start = max(0, j * side_x - side_x // 2)
            x_end = min(width, j * side_x + side_x // 2)
            y_start = max(0, i * side_y - side_y // 2)
            y_end = min(height, i * side_y + side_y // 2)
 
            data_r = image_r[y_start:y_end, x_start:x_end]
            data_g = image_g[y_start:y_end, x_start:x_end]
            data_b = image_b[y_start:y_end, x_start:x_end]
 
            image_point_r[i, j] = np.mean(data_r)
            image_point_g[i, j] = np.mean(data_g)
            image_point_b[i, j] = np.mean(data_b)
 
    # Initialize gain matrices
    rGain = np.zeros((side_num + 1, side_num + 1))
    gGain = np.zeros((side_num + 1, side_num + 1))
    bGain = np.zeros((side_num + 1, side_num + 1))
 
    # Calculate LSC gain
    center_i = side_num // 2
    center_j = side_num // 2
    for i in range(side_num + 1):
        for j in range(side_num + 1):
            rGain[i, j] = image_point_r[center_i, center_j] / image_point_r[i, j]
            gGain[i, j] = image_point_g[center_i, center_j] / image_point_g[i, j]
            bGain[i, j] = image_point_b[center_i, center_j] / image_point_b[i, j]
 
    return rGain, gGain, bGain
 
 
def apply_lsc_correction(image, rGain, gGain, bGain, side_num):
    height, width, _ = image.shape
    side_y = height // side_num
    side_x = width // side_num
 
    corrected_image = np.zeros_like(image, dtype=np.float32)
 
    for i in range(height):
        for j in range(width):
            gainStepX = min(side_num - 1, i // side_y)
            gainStepY = min(side_num - 1, j // side_x)
 
            # Bilinear interpolation
            rGainTmp = (rGain[gainStepX + 1, gainStepY] - rGain[gainStepX, gainStepY]) * (i % side_y) / side_y + \
                       (rGain[gainStepX, gainStepY + 1] - rGain[gainStepX, gainStepY]) * (j % side_x) / side_x + \
                       (rGain[gainStepX + 1, gainStepY + 1] + rGain[gainStepX, gainStepY] - rGain[
                           gainStepX + 1, gainStepY] - rGain[gainStepX, gainStepY + 1]) * \
                       (i % side_y) / side_y * (j % side_x) / side_x + rGain[gainStepX, gainStepY]
 
            gGainTmp = (gGain[gainStepX + 1, gainStepY] - gGain[gainStepX, gainStepY]) * (i % side_y) / side_y + \
                       (gGain[gainStepX, gainStepY + 1] - gGain[gainStepX, gainStepY]) * (j % side_x) / side_x + \
                       (gGain[gainStepX + 1, gainStepY + 1] + gGain[gainStepX, gainStepY] - gGain[
                           gainStepX + 1, gainStepY] - gGain[gainStepX, gainStepY + 1]) * \
                       (i % side_y) / side_y * (j % side_x) / side_x + gGain[gainStepX, gainStepY]
 
            bGainTmp = (bGain[gainStepX + 1, gainStepY] - bGain[gainStepX, gainStepY]) * (i % side_y) / side_y + \
                       (bGain[gainStepX, gainStepY + 1] - bGain[gainStepX, gainStepY]) * (j % side_x) / side_x + \
                       (bGain[gainStepX + 1, gainStepY + 1] + bGain[gainStepX, gainStepY] - bGain[
                           gainStepX + 1, gainStepY] - bGain[gainStepX, gainStepY + 1]) * \
                       (i % side_y) / side_y * (j % side_x) / side_x + bGain[gainStepX, gainStepY]
 
            corrected_image[i, j, 0] = image[i, j, 0] * rGainTmp
            corrected_image[i, j, 1] = image[i, j, 1] * gGainTmp
            corrected_image[i, j, 2] = image[i, j, 2] * bGainTmp
 
    return np.clip(corrected_image, 0, 255).astype(np.uint8)
 
 
if __name__ == "__main__":
    filePath = '../images/LSC.png'
    side_num = 16
    meshON = True
 
    # Load image
    image = cv2.imread(filePath)
 
    # Calculate LSC gain
    rGain, gGain, bGain = calculate_lsc_gain(image, side_num)
 
    # Apply LSC correction
    corrected_image = apply_lsc_correction(image, rGain, gGain, bGain, side_num)
 
    # Display images
    plt.figure(figsize=(12, 6))
    plt.subplot(121)
    plt.imshow(cv2.cvtColor(image, cv2.COLOR_BGR2RGB))
    plt.title('Original Image')
 
    plt.subplot(122)
    plt.imshow(cv2.cvtColor(corrected_image, cv2.COLOR_BGR2RGB))
    plt.title('Corrected Image')
    plt.show()

LSC一般包括两部分：标定和矫正

 

首先问一个问题：为什么需要标定后矫正而不是直接每个点乘与中心点的比值直接进行矫正。

 

有以下几点因素导致：

1、当前采集端分辨率普遍较高，一亿像素的增益表存在芯片里是很占用内存的；

2、在矫正的过程中，像相邻两个像素一个增益1.111和一个增益1.112，重复存下来效率很低，也没必要；

3、芯片里面处理不了浮点数，要存下来的话，就只能提高位深，这样就更占内存了；

 

常见的几种LSC的几种方法如下：

1、radial shading correct

这种方法默认相同半径的点对应的增益是相同的，所以把对应的增益存储在内存中，用到的时候再拿出来，从而完成校正。但是不能把所有的像素的半径都存储起来，所以就通过采样的方式提取特征半径的增益存储到内存中，然后其他半径对应的增益在矫正的时候通过插值算法求出来。

2、mesh shading correct

这种方法是把整幅图像分成n*n个网格，然后针对网格顶点求出矫正的增益，然后把这些顶点的增益储存到内存中，同理其他的点的增益也是通过插值的方法求出

3、多项式拟合

这种方法就是用半径为采样点，然后把这些采样点通过高次拟合成一个高次曲线，然后把高次曲线的参数存储起来，用的时候把半径带入公式就能求出对应的增益值。

同心圆法的流程为：

1. 找到RGB三通道的圆心（一般选择为同一个点）

2. 以同心圆的形状将画面的中心和画面的边缘的三通道乘以不同的增益

具体如下图所示，一般来说考虑shading渐变的曲率从中心到边缘逐渐增大，所以等增益曲线中心稀疏，边缘密集。一般来说lens shading的增益最好不要超过2倍，因为会引入噪声。

下图为mesh shading的矫正方法，同一个方格中的增益一致，mesh的分布也是中心稀疏四角密集。

最后总结一下：同心圆矫正方法的优点是计算量小，缺点是镜头若装配时稍有不对称则矫正失败；网格矫正方法的优点是能够应对各种shanding情况，缺点是运算量大。

 1.2.3. LSC代码

mesh shading correct是当前应用比较广泛的一种LSC，所以这次复现这种算法。在复现算法前，由于这次拿到的raw图是不知道bayer阵列的，所以这里也个大家介绍一种在没有其他工具的帮助下如何确定bayer阵列。

 

基本思路大致为：

• 用待测手机拍一张色卡，同时出.raw和.jpg图
• imshow显示图像时同时会有图像坐标，由此我们大致可以确定红、绿、蓝三个通道的坐标范围
• 在bayer域从奇数行到偶数行，奇数列到偶数列开始输出像素值，色卡红色块中值最大的是R通道，蓝色通达同理，剩下的两个分别是GR和GB。

具体步骤参考如下：

1. 先拍摄一张亮度均匀的raw图，如下图。测下来图片的shading只有30%左右：

2. 将图像分成17 x 13块。（通常sensor出图为 4:3的尺寸，分成17 x 13块可以保证每一块近似于方形）。去除OB后，分通道统计每一块的亮度均值。 

    grid_x = 17
    grid_y = 13
    block_x = int(width / 2 / grid_x)
    block_y = int(high / 2 / grid_y)
    R_otpdata = np.zeros((grid_y, grid_x))
    Gr_otpdata = np.zeros((grid_y, grid_x))
    Gb_otpdata = np.zeros((grid_y, grid_x))
    B_otpdata = np.zeros((grid_y, grid_x))
    data_R = data[::2, ::2]
    data_Gr = data[1::2, ::2]
    data_Gb = data[::2, 1::2]
    data_B = data[1::2, 1::2]
    for i in range(grid_x):
        for j in range(grid_y):
            R_otpdata[j][i] = np.mean(data_R[j * block_y:(j + 1) * block_y, i * block_x:(i + 1) * block_x])
            Gr_otpdata[j][i] = np.mean(data_Gr[j * block_y:(j + 1) * block_y, i * block_x:(i + 1) * block_x])
            Gb_otpdata[j][i] = np.mean(data_Gb[j * block_y:(j + 1) * block_y, i * block_x:(i + 1) * block_x])
            B_otpdata[j][i] = np.mean(data_B[j * block_y:(j + 1) * block_y, i * block_x:(i + 1) * block_x])

如下是Gr通道的亮度统计图（以下的统计图全部以Gr通道为例）：

3. 求出各通道每一块的补偿系数：

    R_otpdata = R_otpdata.max() / R_otpdata
    Gr_otpdata = Gr_otpdata.max() / Gr_otpdata
    Gb_otpdata = Gb_otpdata.max() / Gb_otpdata
    B_otpdata = B_otpdata.max() / B_otpdata

如下是Gr通道需要补偿的系数矩阵：

 4. 将各通道17 x 13的补偿矩阵插值到原本的尺寸，一般来说用cos4次方曲面来拟合插值的效果最好，这里就简单的用双线性插值的方法，最后出来的效果也很好。

    R_otpdata = cv2.resize(R_otpdata,(1296,972),interpolation = cv2.INTER_LINEAR)
    Gr_otpdata = cv2.resize(Gr_otpdata, (1296,972), interpolation=cv2.INTER_LINEAR)
    Gb_otpdata = cv2.resize(Gb_otpdata, (1296,972), interpolation=cv2.INTER_LINEAR)
    B_otpdata = cv2.resize(B_otpdata, (1296,972), interpolation=cv2.INTER_LINEAR)

 如下是插值后的补偿系数矩阵：

5.然后再分通道补偿到原图像中。需要注意不要百分百补偿，原因后面会说到。这里做了一个85%的补偿系数，具体做法感兴趣的同学可以私聊我一起讨论。

补偿后的raw图如下图，shading可以达到85%左右：

1.2.4. 联动方式

Lens shading可不做联动

Color shading可在三种色温下进行差值

1.3. 调试方法

拍摄不同色温的DNP，然后与目标值比较进行校正。

1.4. 测试方法

1.4.1. 客观测试

1.4.1.1. 测试设备

平板光源(DNP)，最好包含三个色温(3000K\5000K\7500K)

如果没有平板光源，拍摄白纸或者白墙也可以，但是测试精度会差。本人实测，即使是DNP均匀度也只有95%-98%，灯墙内壁的均匀度大约为90%-95%。最不济可以自行制作平板光源，需要非常注意随着使用寿命增加光线不均匀的问题。

1.4.1.2. 测试标准

看各自喜好，有的手机厂商要求lens shading在80%以上，行车记录仪或者安防摄像头可适当降低。特别需要注意，如果镜头本身比较差，原始的lens shading为30%，切记不可强行提升过多，因为矫正shading的本质是增加图像边缘部分的增益，这样低照度下会带来非常大的噪声。

有的手机厂商要求color shading的比值在0.9-1.1之间（一般来说R/B的color shading会偏差最大，但不可只测R/B），有的不要求。

1.4.2. 主观测试

1. 拍摄灰墙或者白墙，观察四角发暗的情况

2. 观察图像四角的偏色情况

3. 观察图像四周的跳动噪声

FPGA ISP LensShading

关于LensShading 形成的原因， 参考资料众多，这里重点讲解FPGA上如何实现Mesh LensShading。

首先，给出FPGA实现的系统框图：

在算法设计时，对图像采用分块的方式，每个块中包含一定数目的像素点，实际需要根据处理的分辨率来设计块的大小，这个需要在资源和效果两方面需要考虑：采用较小的块均衡后效果会比较细腻，但是占用大量偏上存储资源保存系数；采用较大的块均衡后效果可能比较粗糙，但是对偏上存储资源要求低。

 

根据输入像素点，分别在行方向和列方向计数，确定当前像素点所在的行块、列块的系数，拼接形成地址，用于索引系数LUT。同时，行计数、列计数在每个块内会有一个局部的坐标，用于索引行权重、列权重。所以FPGA会用到3张查找表，用于计算当前像素点的均衡系数。

如果所示，在对黄色的pixel进行均衡时，知道块索引以后，可以得到当前像素点用于均衡的四个系数（保存于系数LUT），Weight coe1/2/3/4。同时知道黄色pixel在矩阵中相对位置，则可以根据双线性插值原理，计算pixel对应的权重。最后则可以完成对pixel的均衡。

 

理论上mesh的方式要优于radius的方式，sensor毕竟是矩形。另外，海思芯片的ISP也是推荐采用mesh shading的方式。