本文介绍镜头阴影校正(LSC)技术的基本原理及其实现方法。LSC用于解决由镜头光学特性引起的图像边角偏暗问题，包括两种主要校正方法：同心圆法和网格法，并提供了网格法的具体实现代码。
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个网格，然后针对网格顶点求出校正的增益，然后把这些顶点的增益储存到内存中，其它点的增益通过插值的方式求出。
比如图像分成如下的网格：
如下图是每个网格的亮度分布，这里有一个 c o s 4 θ cos^{4}\theta cos4θ的关系。
针对上面的亮度求出的增益图如下：
c o s 4 θ cos^{4}\theta cos4θ的函数如下：
同心圆校正方法的优点是计算量小，缺点是镜头若装配时稍有不对称则校正失败–这个方法可以通过先找到图像的实际中心，然后以实际中心为圆点去校正。
网格法的优点是能够应对各种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*Ny
const 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);
}
}
// 合并raw
general::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));
}
校正后效果：
注1：文中所用图片：
链接1：lens shading原始raw图
链接2：：https://pan.baidu.com/s/1nzdpojJmsgPtKqKzfShqQg?pwd=q45l
提取码：q45l
注2：文中算法有一定缺陷，对网格与网格之间的处理不够细致，放大校正后的图片可以看到有边界感。此算法用来简单理解其原理，更多内容请自行查找其它资料，也欢迎分享相关资料。
扫一扫
7283
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
