ISP之LSC（Lens Shading Correction）

大家好，我是讯享网，很高兴认识大家。

LSC - Lens Shading Correction即镜头暗影校正

一、LSC的意义

众所周知Lens Shading分为Luma Shading和Color Shading。一般来说，物体到Lens中心的距离越远，图像越暗，呈圆形中性对称的方式递减。由于Luma Shading和Color Shading的成因及现象不同，所以我们也分开来分析一下这两种现象。

1、Luma Shading：

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

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

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

Luma Shading的现象如下图：

讯享网

2、Color Shading：

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小。否则将会导致成像照度或色彩问题

Color Shading的图例

二、LSC的实现

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

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

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

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

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

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

1、radial shading correct

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

2、mesh shading correct

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

3、多项式拟合

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

三、代码实现

mesh shading correct是当前应用比较广泛的一种LSC，所以博主这次复现这种算法。

在复现算法前，由于博主这次拿到的raw图是不知道bayer阵列的，所以这里也个大家介绍一种在没有其他工具的帮助下如何确定bayer阵列。

基本思路大致为：

1、用待测手机拍一张色卡，同时出.raw和.jpg图

2、imshow显示图像时同时会有图像坐标，由此我们大致可以确定红、绿、蓝三个通道的坐标范围

3、在bayer域从奇数行到偶数行，奇数列到偶数列开始输出像素值，色卡红色块中值最大的是R通道，蓝色通达同理，剩下的两个分别是GR和GB。

 imgrgb = cv2.imread('/Users/ISP/LSC/color/Camera/IMG__.jpg') #bgr转为rgb imgrgb = imgrgb[:, :, [2, 1, 0]] plt.imshow(imgrgb) plt.show() # 红色块坐标范围(1600, 1777)-(1753, 1939) # 绿色块坐标范围（1600，1516）- （1753，1678） # 蓝色块坐标范围（1600，1255）-（1753，1417） img_color = np.fromfile(r'/Users/ISP/LSC/color/Camera/IMG__.dng',dtype = 'uint16') L = len(img_color) a = 4032 * 3024 img_color16 = img_color[L-a:L].reshape(3024,4032) #矩阵奇数行到偶数行，奇数列到偶数列 #python从0开始所以从偶数行到奇数行，从偶数列到奇数列 #按照计算结果bayer阵列为BGGR #红块[132,211,194,258] print(img_color16[1672:1674, 1856:1858]) #蓝[[228,173,165,94] print(img_color16[1672:1674, 1356:1358]) #np.savetxt('/Users/ISP/LSC/color/Camera/new.csv', img_color8, delimiter = ',')

讯享网

确认完bayer阵列，开始做LSC

讯享网 #拆分通道，我切割的比较大，大家在做切割的时候选择高/宽整除的数即可 side_num = 8 img1 = np.fromfile(r'/Users/Desktop/ISP/LSC/IMG__.dng',dtype = 'uint16') a = len(img1) b = 4032 * 3024 #dng与raw的差别是，dng有头部信息，需要剪掉头部信息才是真实像素值，原图为10bit图像 img = img1[a-b:a].reshape(3024,4032) B = img[::2, ::2] GB = img[::2, 1::2] GR = img[1::2, ::2] R = img[1::2, 1::2] R_max = np.max(R) GR_max = np.max(GR) GB_max = np.max(GB) B_max = np.max(B) Hblocks = int(R.shape[0]/side_num) Wblocks = int(R.shape[1]/side_num) R_LSC_data = np.zeros((Hblocks, Wblocks)) B_LSC_data = np.zeros((Hblocks, Wblocks)) GR_LSC_data = np.zeros((Hblocks, Wblocks)) GB_LSC_data = np.zeros((Hblocks, Wblocks)) Rgain_tmp = np.zeros((Hblocks, Wblocks)) Grgain_tmp = np.zeros((Hblocks, Wblocks)) Gbgain_tmp = np.zeros((Hblocks, Wblocks)) Bgain_tmp = np.zeros((Hblocks, Wblocks)) img_new = np.zeros((3024, 4032)) gain = np.zeros((3024, 4032)) #标定 for y in range(0, R.shape[0], side_num): for x in range(0, R.shape[1], side_num): block_y_num = int(y / side_num) block_x_num = int(x / side_num) R_LSC_data[block_y_num, block_x_num] = R[y:y + side_num, x: x + side_num].mean() GR_LSC_data[block_y_num, block_x_num] = GR[y:y + side_num, x: x + side_num].mean() GB_LSC_data[block_y_num, block_x_num] = GB[y:y + side_num, x: x + side_num].mean() B_LSC_data[block_y_num, block_x_num] = B[y:y + side_num, x: x + side_num].mean() Rgain_tmp = R_max/R_LSC_data Grgain_tmp = GR_max/GR_LSC_data Gbgain_tmp = GB_max/GB_LSC_data Bgain_tmp = B_max/B_LSC_data Rgain = np.zeros((R.shape[0],R.shape[1])) GRgain = np.zeros((R.shape[0],R.shape[1])) GBgain = np.zeros((R.shape[0],R.shape[1])) Bgain = np.zeros((R.shape[0],R.shape[1])) #插值矫正，调用opencv简单实现 #Rgain = cv2.resize(Rgain_tmp,(R.shape[1],R.shape[0]) ,interpolation=cv2.INTER_CUBIC) #GRgain = cv2.resize(Grgain_tmp, (R.shape[1],R.shape[0]),interpolation=cv2.INTER_CUBIC) #GBgain = cv2.resize(Gbgain_tmp, (R.shape[1],R.shape[0]),interpolation=cv2.INTER_CUBIC) #Bgain = cv2.resize(Bgain_tmp, (R.shape[1],R.shape[0]),interpolation=cv2.INTER_CUBIC) #ISP芯片中没有opencv，用手动插值实现，我们本次尝试用双线性插值 for j in range(R.shape[1]): for i in range(R.shape[0]): src_x = i / side_num src_y = j / side_num x = int(math.floor(src_x)) y = int(math.floor(src_y)) u = src_x - x v = src_y - y if x == 188: x = 187 if y == 251: y = 250 Rgain[i, j] = (1-u) * (1-v) * Rgain_tmp[x,y] + (1-u) * v*Rgain_tmp[x,y + 1] + u*(1-v)*Rgain_tmp[x + 1,y] +\ u*v*Rgain_tmp[x + 1,y + 1] GRgain[i, j] = (1 - u) * (1 - v) * Grgain_tmp[x, y] + (1 - u) * v * Grgain_tmp[x, y + 1] + u * (1 - v) * \ Grgain_tmp[x + 1, y] + u * v * Grgain_tmp[x + 1, y + 1] GRgain[i, j] = (1 - u) * (1 - v) * Gbgain_tmp[x, y] + (1 - u) * v * Gbgain_tmp[x, y + 1] + u * (1 - v) * \ Gbgain_tmp[x + 1, y] + u * v * Gbgain_tmp[x + 1, y + 1] GBgain[i, j] = (1 - u) * (1 - v) * Grgain_tmp[x, y] + (1 - u) * v * Grgain_tmp[x, y + 1] + u * (1 - v) * \ Grgain_tmp[x + 1, y] + u * v * Grgain_tmp[x + 1, y + 1] Bgain[i, j] = (1 - u) * (1 - v) * Bgain_tmp[x, y] + (1 - u) * v * Bgain_tmp[x, y + 1] + u * (1 - v) * \ Bgain_tmp[x + 1, y] + u * v * Bgain_tmp[x + 1, y + 1] gain[::2, ::2] = Rgain gain[::2, 1::2] = GRgain gain[1::2, ::2] = GBgain gain[1::2, 1::2] = Bgain img_new = img*gain img_new = img_new/4 plt.subplot(1, 2, 1) plt.imshow(img, 'gray') plt.subplot(1, 2, 2) plt.imshow(img_new, 'gray') plt.show() cv2.imwrite('/Users/ISP/LSC/LSC_correct.jpg',img_new)

四、结果

参考文献：

关于CRA可以参考这片文章深入浅出Sensor CRA(主光角)|极客笔记