ISP-LSC 镜头阴影校正

参考：

https://zhuanlan.zhihu.com/p/389334269
https://blog.csdn.net/xiaoyouck/article/details/77206505
https://www.cnblogs.com/wnwin/p/11805901.html
http://kb.colorspace.com.cn/kb/2022/09/05/isp-%E9%95%9C%E5%A4%B4%E9%98%B4%E5%BD%B1%E6%A0%A1%E6%AD%A3%EF%BC%88lsc%EF%BC%89/

LSC(Lens Shading Correction)镜头阴影校正。

Shading细分为Lens Shading(Luma Shaing)和Color Shading(chrom Shading)。

1. Shading产生原因及影响

由于Lens的光学特性，sensor影像区的边缘区域接收的光强比中心小，造成中心和四角亮度不一致的现象，并且镜头本身是一个凸透镜，由于凸透镜的原理，中心的感光必然比周边多。

如下图，绿色和蓝色的光强在进入镜头前是一致的，由于Lens的特性，在边缘的绿色光线会有一部分被遮挡，sensor边缘部分能捕捉到的光信号就比较少；

另外由于绿色部分光线经过的距离比较长，光的衰减也要比蓝色部分的衰减大，也导致了到达边缘部分的光信号的强度弱；

以上来自：https://zhuanlan.zhihu.com/p/389334269

在这里插入图片描述

Lens Shading会造成图像边角偏暗，即暗角。

在这里插入图片描述

color Shading

各种颜色的波长不同，经过透镜折射后，折射的角度也不一样，就会造成color Shading的现象，另外由于CRA的原因也会导致shading现象，如下：

在这里插入图片描述

Color Shading 中心和四周颜色不一致，表现出来一般是中心或者四周偏色。

在这里插入图片描述

2. Shading校正-LSC

在ISP Pipeline中，Shading一般在OB和DPC后面。另外需要注意，如果3A的统计数据是在shading校正之后获取的，那么shading校正结果会影响3A的统计数据。

2.1 校正方法

LSC校正目前主流的方法有两种：同心圆法和网格法。

同心圆法：

找到RGB三通道的圆心；
以同心圆的形状将画面的中心和画面的边缘的三通道乘以不同的增益。

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

在这里插入图片描述

网格法：

也叫做mesh shading correction，把整幅图像分成m*n个网格，然后针对网格顶点求出校正的增益，然后把这些顶点的增益储存到内存中，其它点的增益通过插值的方式求出。

比如图像分成如下的网格：

在这里插入图片描述

如下图是每个网格的亮度分布，这里有一个 $cos^{4}\theta$ 的关系。

在这里插入图片描述

针对上面的亮度求出的增益图如下：

在这里插入图片描述

$cos^{4}\theta$ 的函数如下：

同心圆校正方法的优点是计算量小，缺点是镜头若装配时稍有不对称则校正失败–这个方法可以通过先找到图像的实际中心，然后以实际中心为圆点去校正。

网格法的优点是能够应对各种shading情况，缺点是运算量大。

网格法校正代码流程：

整张raw图像分为四通道；
对四通道图像划分为m*n个网格，并求出每个网格的均值；
根据每个网格均值求出该网格对应的增益；
根据得到的增益再利用插值算法得到每个像素校正后的值；

template <typename T>
bool lsc(const T *src, T *dst, int width, int height)
{// 拆分为四通道int sub_width = width / 2;int sub_height = height / 2;T *pr, *pb, *pgr, *pgb;pr = new T[sub_width * sub_height];pb = new T[sub_width * sub_height];pgr = new T[sub_width * sub_height];pgb = new T[sub_width * sub_height];if (!pr || !pb || !pgr || !pgb){printf("allocate channel buf failed!\n");return false;}memset(pr, 0, sizeof(pr[0]) * sub_width * sub_height);memset(pb, 0, sizeof(pb[0]) * sub_width * sub_height);memset(pgr, 0, sizeof(pgr[0]) * sub_width * sub_height);memset(pgb, 0, sizeof(pgb[0]) * sub_width * sub_height);general::raw_split_to_4channel(src, pr, pgr, pgb, pb, width, height);// 整张图像划分为17*13 block，并获取平均值int grid_width, grid_height; //每个网格的像素个数const int Nx = 17; // 网格个数 Nx*Nyconst int Ny = 13;const int nx = Nx + 1;const int ny = Ny + 1;grid_width = floor(sub_width * 1.0 / (Nx));grid_height = floor(sub_height * 1.0 / (Ny));float rave[Ny + 1][Nx + 1],grave[Ny + 1][Nx + 1],gbave[Ny + 1][Nx + 1],bave[Ny + 1][Nx + 1];memset(rave, 0, sizeof(rave));memset(grave, 0, sizeof(grave));memset(gbave, 0, sizeof(gbave));memset(bave, 0, sizeof(bave));int sx, sy, ex, ey;sx = sy = ex = ey = 0;for (int i = 0; i <= Ny; i++){for (int j = 0; j <= Nx; j++){sx = j * grid_width - grid_width / 2;ex = j * grid_width + grid_width / 2;sy = i * grid_height - grid_height / 2;ey = i * grid_height + grid_height / 2;if (i == Ny && ey != sub_height)ey = sub_height;if (j == Nx && ex != sub_width)ex = sub_width;sx = sx < 0 ? 0 : sx;sy = sy < 0 ? 0 : sy;rave[i][j] = general::get_average_roi(pr, sub_width, sx, sy, ex, ey);grave[i][j] = general::get_average_roi(pgr, sub_width, sx, sy, ex, ey);gbave[i][j] = general::get_average_roi(pgb, sub_width, sx, sy, ex, ey);bave[i][j] = general::get_average_roi(pb, sub_width, sx, sy, ex, ey);}}// 获取每个通道的均值最大值float max[4] = {0, 0, 0, 0};for (int i = 0; i <= Ny; i++){for (int j = 0; j <= Nx; j++){max[0] = (max[0] < rave[i][j] ? rave[i][j] : max[0]);max[1] = (max[1] < grave[i][j] ? grave[i][j] : max[1]);max[2] = (max[2] < gbave[i][j] ? gbave[i][j] : max[2]);max[3] = (max[3] < bave[i][j] ? bave[i][j] : max[3]);}}// 计算每个通道的增益float rgain[ny][nx], grgain[ny][nx],gbgain[ny][nx], bgain[ny][nx];memset(rgain, 0, sizeof(rgain));memset(grgain, 0, sizeof(grgain));memset(gbgain, 0, sizeof(gbgain));memset(bgain, 0, sizeof(bgain));for (int i = 0; i <= Ny; i++){for (int j = 0; j <= Nx; j++){rgain[i][j] = max[0] / float(rave[i][j]);grgain[i][j] = max[1] / float(grave[i][j]);gbgain[i][j] = max[2] / float(gbave[i][j]);bgain[i][j] = max[3] / float(bave[i][j]);}}// 计算最终值float gaintmp = 0;int gainx, gainy;gainx = gainy = 0;int tmp_grid_width, tmp_grid_height;int tmp_x, tmp_y;float tmp = 0;int curgain = 0;T *prdst, *pbdst, *pgrdst, *pgbdst;prdst = new T[sub_width * sub_height];pbdst = new T[sub_width * sub_height];pgrdst = new T[sub_width * sub_height];pgbdst = new T[sub_width * sub_height];if (!prdst || !pbdst || !pgrdst || !pgbdst){printf("allocate channel buf failed!\n");return false;}memset(prdst, 0, sizeof(prdst[0]) * sub_width * sub_height);memset(pbdst, 0, sizeof(pbdst[0]) * sub_width * sub_height);memset(pgrdst, 0, sizeof(pgrdst[0]) * sub_width * sub_height);memset(pgbdst, 0, sizeof(pgbdst[0]) * sub_width * sub_height);tmp = 0;for (int y = 0; y < sub_height; y++){for (int x = 0; x < sub_width; x++){gainy = floor(float(y) / grid_height);gainy = (gainy > Ny - 1) ? (Ny - 1) : gainy;gainx = floor(float(x) / grid_width);gainx = (gainx > Nx - 1) ? (Nx - 1) : gainx;// 			if (x == 2103 && y == 1557)// 				tmp = 1;// 插值得到每个像素校正后的值// f(x,y) = [f(1,0)-f(0,0)]*x +// [f(0,1) - f(0,0)]*y +// [f(1,1)+f(0,0)-f(0,1)-f(1,0)]*xy +// f(0,0)gaintmp = (rgain[gainy][gainx + 1] - rgain[gainy][gainx]) * (x - gainx * grid_width) / grid_width +(rgain[gainy + 1][gainx] - rgain[gainy][gainx]) * (y - gainy * grid_height) / grid_height +(rgain[gainy + 1][gainx + 1] + rgain[gainy][gainx] - rgain[gainy + 1][gainx] - rgain[gainy][gainx + 1]) * (x - gainx * grid_width) / grid_width * (y - gainy * grid_height) / grid_height +rgain[gainy][gainx];prdst[x + y * sub_width] = T(float(pr[x + sub_width * y]) * gaintmp);gaintmp = (grgain[gainy][gainx + 1] - grgain[gainy][gainx]) * (x - gainx * grid_width) / grid_width +(grgain[gainy + 1][gainx] - grgain[gainy][gainx]) * (y - gainy * grid_height) / grid_height +(grgain[gainy + 1][gainx + 1] + grgain[gainy][gainx] - grgain[gainy + 1][gainx] - grgain[gainy][gainx + 1]) * ((x - gainx * grid_width) / grid_width) * ((y - gainy * grid_height) / grid_height) +grgain[gainy][gainx];pgrdst[x + y * sub_width] = T(float(pgr[x + sub_width * y]) * gaintmp);gaintmp = (gbgain[gainy][gainx + 1] - gbgain[gainy][gainx]) * (x - gainx * grid_width) / grid_width +(gbgain[gainy + 1][gainx] - gbgain[gainy][gainx]) * (y - gainy * grid_height) / grid_height +(gbgain[gainy + 1][gainx + 1] + gbgain[gainy][gainx] - gbgain[gainy + 1][gainx] - gbgain[gainy][gainx + 1]) * ((x - gainx * grid_width) / grid_width) * ((y - gainy * grid_height) / grid_height) +gbgain[gainy][gainx];pgbdst[x + y * sub_width] = T(float(pgb[x + sub_width * y]) * gaintmp);gaintmp = (bgain[gainy][gainx + 1] - bgain[gainy][gainx]) * (x - gainx * grid_width) / grid_width +(bgain[gainy + 1][gainx] - bgain[gainy][gainx]) * (y - gainy * grid_height) / grid_height +(bgain[gainy + 1][gainx + 1] + bgain[gainy][gainx] - bgain[gainy + 1][gainx] - bgain[gainy][gainx + 1]) * ((x - gainx * grid_width) / grid_width) * ((y - gainy * grid_height) / grid_height) +bgain[gainy][gainx];pbdst[x + y * sub_width] = T(float(pb[x + sub_width * y]) * gaintmp);}}// 合并rawgeneral::channels4_to_raw(dst, prdst, pgrdst, pgbdst, pbdst, width, height);return true;
}/// @brief raw图拆分为四通道图/// @tparam T 传入的src及p1/p2/p3/p4 buf的类型/// @param src in 原始raw buf，pixel raw格式/// @param p1p2p3p4 out 拆分后得到的四通道的buf，不分通道，位置如下：///                 |p1|p2|///                  —— ——///                 |p3|p4|/// @param width in raw宽/// @param height in raw高template <typename T>void raw_split_to_4channel(const T *src, T *p1, T *p2, T *p3, T *p4, int width, int height){int index = 0;for (int y = 0; y < height; y += 2){for (int x = 0; x < width; x += 2){p1[index] = src[x + y * width];p2[index] = src[x + y * width + 1];p3[index] = src[x + (y + 1) * width];p4[index] = src[x + (y + 1) * width + 1];index++;}}}/// @brief 四通道图合并为原始raw/// @tparam T 传入的src及p1/p2/p3/p4 buf的类型/// @param dst out 合并后的完整raw图/// @param p1p2p3p4 in 四通道的buf，不分通道，位置如下：///                 |p1|p2|///                  —— ——///                 |p3|p4|/// @param width in raw宽/// @param height in raw高template <typename T>void channels4_to_raw(T *dst, T *p1, T *p2, T *p3, T *p4, int width, int height){int index = 0;for (int y = 0; y < height; y += 2){for (int x = 0; x < width; x += 2){dst[x + y * width] = p1[index];dst[x + y * width + 1] = p2[index];dst[x + (y + 1) * width] = p3[index];dst[x + (y + 1) * width + 1] = p4[index];index++;}}}/// @brief 获取raw buf中某个roi区域的平均值/// @tparam T 传入的src及p1/p2/p3/p4 buf的类型/// @param src in raw buf/// @param w in raw 宽/// @param sx in roi的位置，start x/// @param sy in roi的位置，start y/// @param ex in roi的位置，end x/// @param ey in roi的位置，end y/// @return raw buf中该roi的均值template <typename T>double get_average_roi(const T *src, int w, int sx, int sy, int ex, int ey){double ave = 0.0;if (sx == ex || sy == ey)return ave;ave = 0;for (int y = sy; y < ey; ++y){for (int x = sx; x < ex; ++x){ave += src[x + y * w];}}return ave / ((ex - sx) * (ey - sy));}