RAW域处理算法之LSC

Posted 2023-06-22 lanlancky

RAW域处理算法之LSC

实际应用中，由于具体场景的需要以及成本的考虑，摄像机会搭配不同镜头。镜头校正是指针对由于镜头原因引入的成像误差进行的校正。

镜头阴影校正（Lens shading correction，LSC）

由于镜头/微镜头的光学构造，导致了经过镜头/微镜头进入 sensor 的光线中间多四周少，进而导致了数字图像出;现中间区域亮，四周区域暗(shading)的情况。通过暗角校正期望校正后的图像的四周和中心亮度一致。
校正算法原理：测定shading 的分布情况，采用多项式拟合或同心圆补偿的方法对其校正。

镜头阴影有两种表现形式，分别是

• Luma shading，又称vignetting，指由于镜头通光量从中心向边缘逐渐衰减导致画面边缘亮度变暗的现象。如下图左所示。
• Chroma shading，指由于镜头对不同波长的光线折射率不同引起焦平面位置分离导致图像出现伪彩的现象。如下图右所示。

Luma shading

• 画面边缘镜头能量衰减：如图所示，蓝色和绿色用相同的数量线条表示能量，中心位置的蓝色几乎所有能量都能达到最右侧的的成像单元，但是边缘的绿色由于有一定角度射入，经过镜头的折射，有一部分光（最上方的几条绿色线条）就没法达到成像单元，因此成像单元中心的能量就会比边缘的大，变现在亮度上就是亮度向边缘衰减变暗。通常镜头的衰减符合。
• 所以Luma shading的主要原因是镜头中心到边缘的能量衰减导致的，由于镜头中都会存在多处光阑，当入射光线偏离光轴角度较大时，部分光线就会被光阑遮挡而不能参与成像，因此越靠近sensor边缘的像素接收到的曝光量就越低。

• 边缘像素微透镜和感光面的错位：这个问题在手机sensor 上通常会更严重一些，因此设计手机sensor 的厂家会采取一些特定的方法去缓解这个问题。一种常用的方法是在微透镜上做文章，即从中心像素开始，微透镜的尺寸略小于感光面的面积一点点，这样越往边缘微透镜与感光面之间的错位就越大，刚好可以补偿入射光线角度增大导致的焦点偏移，使光线可以更好地聚焦到感光面上，如下图所示。

焦点错位导致能量损失

不过深入研究会发现，这个补偿办法其实也是有局限的，如果sensor采用的是下图左所示的FSI工艺（前照式），从像素微观结构来看，当入射光线角度比较大时，会有较多光线与像素中的金属布线层发生吸收、散射从而产生损失，单纯移动微透镜的位置并不能有效解决这个问题。但是，如果像素采用的是右图所示的BSI工艺（背照式），因为布线层在硅片的另外一侧，所以光线损失会少，补偿效果更加有效。

Vignetting是由镜头引起的现象，所以LSC校正也是针对特定镜头的。如果产品的适配镜头发生变化，原则上需要重新进行LSC校正。另外，Vignetting现象在sensor靶面较大、镜头焦距较短时表现更加明显。采用非球面镜头通常可以改善vignetting。

Chroma shading 原理

镜头对不同波长的光线折射率不同会导致色差问题，即不同波长的焦点在空间上不重合，导致焦平面分裂为三个不完全重合的曲面，这会破坏图像的白平衡，使图像出现伪彩，如下图所示。根据sensor所处的前后位置不同，伪彩可能偏红也可能偏蓝。

镜头带来的color shading主要是因为不同颜色的光的折射率不同，导致白光经过镜头后达到成像面时不同颜色的光的位置不同导致偏色。当然偏色还会和CRA有关，但是一般镜头选型的时候都会注意和sensor的CRA进行匹配，Chroma shading 一般主要通过镜头选型来控制其影响。

如果ISP支持chroma shading校正，则可以通过标定三个颜色平面的增益来修正。为了控制标定表格的存储空间，通常只标定MxN个关键点，任意位置处的像素增益可以使用相邻四个标定关键点通过双线性插值的方法动态计算得到。在一个典型的ISP实现中，可以取M=N=17即可得到令人满意的效果。

LSC矫正方法

LSC的本值就是能量有衰减，反过来为了矫正就用该点的像素值乘以一个gain值，让其恢复到衰减前的状态，所以矫正的本质就是找到这个gain值。

从目前的矫正方法来看个人觉得可以分成三大类，但目前方法1和方法2是使用最多的。

1. 储存增益法
2. 多项式拟合法
3. 自动矫正法

储存增益法

“cos4定律”说明，减少的光线与增大角度余弦值的四次方是成比例关系的。另外，在某些镜头设计中，镜头可能本身就会阻挡一部分光线(称为“晕光”)，这也会引起阴影现象。所以，即使微镜头设计可以最小化短镜头的阴影现象，此种现象还是会多多少少地存在。因此后端ISP一般都会有LSC功能。

上面有提到衰减符合cos(θ)的四次方规律，而θ在三维空间对各个方向是一致的，所以各个方向的衰减如下图

图中相同颜色可以理解成亮度是一样的，也就是图中红色一圈圈的像素需要的增益是一样的，所以就可以用半径为变量来求出不同半径像素需要的增益。然后把半径对应的增益值储存在内存中，到了要用的时候再拿出来用，从而完成矫正。但是不可能把所有像素的半径都存储起来，所以就通过采样的方式提取特征半径的增益存储到内存，然后其他半径对应的增益在矫正的时候通过插值算法求出来。这种方式对内存的硬件要求就低了。这就是radial shading correct。

mesh shading correct

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

如图是上图分成网格后，每个网格亮度的分布，可以看出和cos(θ)的四次方很接近，然后针对这样的网格亮度求出增益如下图，刚好和亮度分布相反

多项式拟合

多项式拟合的方式就是用半径为采样点，然后把这些采样点通过高次拟合的方式拟合成一个高次曲线，然后把高次曲线的参数储存起来，用的时候把半径带入公式就能求出对应的gain值用于矫正。基于FPGA的镜头光照度补充算法研究及实现中提到可以采用拉格朗日插值法进行插值，具体公式如下图

其他通道也能这么拟合出曲线从而完成各个通道的矫正。

针对光学中心可能不是图像中心的问题，有的论文就提出先在图像中找出光学中心，然后以光学中心为真实中心完成标定，最后也是通过光学中心来求半径带入公式求gain值。

LSC算法具体实现

标定代码如下

clc;clear;close all;

% --------parameters of calibretion------------
filePath = \'../images/lscRefImg.jpg\';
side_num = 17;
meshON = 1;
% ---------------------------------------------

image = imread(filePath);
[height, width] = size(image);
side_y = floor(height/side_num);
side_x = floor(width/side_num);
h = imshow(image);
if meshON
    for i = 1:side_num-1
        line([i*side_x, i*side_x], [1, height], \'color\', \'r\');
        line([1, width], [i*side_y, i*side_y], \'color\', \'r\');
    end        
end
title(\'refImg\');

%% compress resolution
image_point = zeros(side_num+1,side_num+1);
for i = 0:side_num
    for j = 0:side_num
        x_clip = floor([j*side_x - side_x/2, j*side_x + side_x/2]);
        y_clip = floor([i*side_y - side_y/2, i*side_y + side_y/2]);
        % make sure that the last point on the edge
        if(i==side_num && y_clip(2) ~= height) 
            y_clip(2) = height;
        end
        if(j==side_num && x_clip(2) ~= width) 
            x_clip(2) = width;
        end
        x_clip(x_clip<1) = 1;
        x_clip(x_clip>width) = width;
        y_clip(y_clip<1) = 1;
        y_clip(y_clip>height) = height;
        data_in = image(y_clip(1):y_clip(2), x_clip(1):x_clip(2));
        image_point(i+1,j+1) = mean(mean(data_in));
    end
end

Gain = zeros(side_num+1,side_num+1);

%% caculate lsc luma gain
for i = 1:side_num+1
    for j = 1:side_num+1
        Gain(i,j) = image_point(uint8(side_num/2) +1, uint8(side_num/2) +1) / image_point(i,j);
    end
end
save(\'./data/Gain.mat\', \'Gain\');

clc, clear, close all;
% --------parameters of correction------------
filePath = \'../images/lscRefImg.jpg\';
side_num = 17;
% --------------------------------------------

% --------load data---------------------------
% load org image
image = imread(filePath);
[height, width] = size(image);
sideX = floor(height/side_num);
sideY = floor(width/side_num);

% load gain of each channel
load(\'./data/Gain.mat\');

% --------------correction-------------------
disImg = zeros(size(image));
gainStepX = 0;
gainStepY = 0;
gainTab = zeros(size(image));
for i = 1:height
    for j = 1:width
        gainStepX = floor(i / sideX) + 1;
        if gainStepX > 16
            gainStepX = 16;
        end
        gainStepY = floor(j / sideY) + 1;
        if gainStepY > 16
            gainStepY = 16;
        end
        % get tht gain of the point by interpolation(Bilinear interpolation)
        % f(x,y) = [f(1,0)-f(0,0)]*x+[f(0,1)-f(0,0)]*y+[f(1,1)+f(0,0)-f(1,0)-f(0,1)]*xy+f(0,0)
        gainTab(i, j) = (Gain(gainStepX+1, gainStepY) - Gain(gainStepX, gainStepY)) 
　　　　　　　　　　　　　　　　* (i - (gainStepX - 1) * sideX)/sideX +...
                        (Gain(gainStepX, gainStepY+1) - Gain(gainStepX, gainStepY)) 
　　　　　　　　　　　　　　　　* (j - (gainStepY - 1) * sideY)/sideY +...
                        (Gain(gainStepX+1, gainStepY+1) + Gain(gainStepX, gainStepY) 
　　　　　　　　　　　　　　　　- Gain(gainStepX+1, gainStepY)- Gain(gainStepX, gainStepY + 1)) *...
                        (i - (gainStepX - 1) * sideX)/sideX * (j - (gainStepY - 1) * sideY)/sideY + ...
 　　　　　　　　　　　　　　Gain(gainStepX, gainStepY);
    end
end
disImg = double(image) .* gainTab;
figure();
subplot(121);imshow(image);title(\'org image\');
subplot(122);imshow(uint8(disImg));title(\'corrected image\');