isp算法学习-LSC
700 Words|Read in about 4 Min|本文总阅读量次
本文是isp算法学习的第三篇章，主要是讲述LSC原理和算法。
1产生原因
shading分成Luma shading和Color shading两种。
1.1Luma shading
影响亮度shading的主要跟镜头有关。用通俗易懂的方式解释，以左边的图为例，蓝色和绿色有相同的能量的光，通过一些列的透镜，蓝色的中心光能量会聚焦到中心，且能量衰减较小，而边缘有一定角度绿色的光会聚焦到边缘，且能量衰减较大。
通常能量衰减公式为
$$ cos^{4}(\theta) $$
可以看到中间零度的时候，能量几乎没有衰减；四周角度大的时候，根据上述公式，发现能量会比中间衰减来的大。
1.2Color shading
光线过来，经过一个棱镜，不同颜色的光，折射率的不同，光会发生色散现象。
经过透镜后有一个平面的话，在平面上就会看到不同颜色的分离，由于颜色的分离，导致中间和四周的颜色会不一致。中心颜色会偏红，四周颜色会偏暗绿。
2校正方法
2.1LUT（显示查找表）
因为镜头和光线折射率的原因，中心点的值为255，那么四周的值可能就只有200，那么需要用一种算法将四周的值乘以一个增益值，是所有的画面像素点的值都能够归一化。通常会取中心像素值最大的点（或者取中心区域的一些点取平均）为目标值。
如果一张图片的像素是800w，如果使用每一个点那么对于芯片整体的内存或者是性能要求会特别高。这个时候我们采取将整个图形分成若干块，存储若干块的顶点即可，这种方式就是查找表方法。
查找表的方式有两种网格法和半径法。
2.1.1网格法
我们对整个图形采用网格法来替代所有像素点的方法，只需要存储网格上的顶点即可。举例以17*17的网格，来存储对应的顶点增益数据，对于800w的数据来说，会大大减轻计算量。
如果说需要求出每一个网格内某一个像素对应的增益值，采用插值的方法来求其增益值即可，不需要再逐个点求增益值。
2.1.2半径法
根据上面公式得出，每个方向都符合cos4(θ)。
相同半径，θ不变，这半径所在的一圈增益也都是相同的。
我取一条线上不同半径的点，每个不同半径上的点保存到对应的表格中。
如果不同半径对应一种增益值，那么任意再图上的点到中心举例 $$ d=\sqrt{(x - x_{c})^2 + (y - y_{c})^2} $$ 但是这种方法实际中很少用到，比较理想化。不能保证镜头安装，工艺生产等方面是否会出现偏差，导致镜头中心并非是画面中心，那么校正会越来越错。
2.2多项式拟合法
2.2.1以图像中心为标定中心
通过拉格朗日插值法，进行高次曲线拟合的增益曲线。
2.2.2不以图像中心为标定中心
镜头安装，工艺生产等方面是否会出现偏差，图像中心并不是我们认为的最亮的点，所以不能简单的以图像中心为标定点。
利用迭代的方式，先垂直的找然后水平的找。定义一个范围TH，在这个范围里面定义一个gain1，超过这个范围的定义gain2，得到两个二项式的拟合曲线，尽可能把真实的gain曲线包住。
3方法实现
可以看到，图像左下角是对应的亮度曲线图，右下角是对应的gain值曲线图，何如对应的gain值后，LSC后的图像就为右上角所示。
3.1标定
均匀光下，采集一张灰色（18%的灰）图像或者滤光片即可。
每个像素都除以整个像素中最大值的80%作为增益
1function LSCCalibrationM(path)
2lscRefImg = double(imread(path));
3tmp = ones(size(lscRefImg));
4corTab = (tmp./lscRefImg) * 0.8 * max(max(lscRefImg));
5save('src/corTab.mat', 'corTab');
6end
3.2直接校正
1%% --------------------------------
2lscRefImg = double(imread('images/lscRefImg.jpg'));
3load('src/data/corTab.mat')
4corImg = uint8(lscRefImg .* corTab);
5figure;
6subplot(121);imshow(uint8(lscRefImg));title('org');
7subplot(122);imshow(corImg);title('corrected');
3.3划分网格校正
1% --------parameters of correction------------
2filePath = 'images/lscRefImg.jpg';
3side_num = 17;
4% --------------------------------------------
5
6% --------load data---------------------------
7% load org image
8image = imread(filePath);
9[height, width] = size(image);
10sideX = floor(height/side_num);
11sideY = floor(width/side_num);
12
13% load gain of each channel
14load('./src/data/Gain.mat');
15
16% --------------correction-------------------
17disImg = zeros(size(image));
18gainStepX = 0;
19gainStepY = 0;
20gainTab = zeros(size(image));
21for i = 1:height
22 for j = 1:width
23 gainStepX = floor(i / sideX) + 1;
24 if gainStepX > 16
25 gainStepX = 16;
26 end
27 gainStepY = floor(j / sideY) + 1;
28 if gainStepY > 16
29 gainStepY = 16;
30 end
31 % get tht gain of the point by interpolation(Bilinear interpolation)
32 % 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)
33 gainTab(i, j) = (Gain(gainStepX+1, gainStepY) - Gain(gainStepX, gainStepY)) * (i - (gainStepX - 1) * sideX)/sideX +...
34 (Gain(gainStepX, gainStepY+1) - Gain(gainStepX, gainStepY)) * (j - (gainStepY - 1) * sideY)/sideY +...
35 (Gain(gainStepX+1, gainStepY+1) + Gain(gainStepX, gainStepY) - Gain(gainStepX+1, gainStepY)- Gain(gainStepX, gainStepY + 1)) *...
36 (i - (gainStepX - 1) * sideX)/sideX * (j - (gainStepY - 1) * sideY)/sideY + Gain(gainStepX, gainStepY);
37 end
38end
39disImg = double(image) .* gainTab;
40
41figure();
42subplot(121);imshow(image);title('org image');
43subplot(122);imshow(uint8(disImg));title('corrected image');
3.4RGB三通道划分
正式版本应该是bayer格式，四通道，改成四通道即可。
1% --------parameters of correction------------
2filePath = 'images/lsc.bmp';
3side_num = 16;
4% --------------------------------------------
5
6% --------load data---------------------------
7% load org image
8image = imread(filePath);
9[height, width, channel] = size(image);
10sideX = floor(height/side_num);
11sideY = floor(width/side_num);
12
13image_r = image(:,:,1);
14image_g = image(:,:,2);
15image_b = image(:,:,3);
16
17% load gain of each channel
18load('./src/data/rGain.mat');
19load('./src/data/gGain.mat');
20load('./src/data/bGain.mat');
21% --------------correction-------------------
22disImg = zeros(size(image));
23gainStepX = 0;
24gainStepY = 0;
25for i = 1:height
26 for j = 1:width
27 gainStepX = floor(i / sideX) + 1;
28 if gainStepX > 16
29 gainStepX = 16;
30 end
31 gainStepY = floor(j / sideY) + 1;
32 if gainStepY > 16
33 gainStepY = 16;
34 end
35 % get tht gain of the point by interpolation(Bilinear interpolation)
36 % 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)
37 rGainTmp = (rGain(gainStepX+1, gainStepY) - rGain(gainStepX, gainStepY)) * (i - (gainStepX - 1) * sideX)/sideX +...
38 (rGain(gainStepX, gainStepY+1) - rGain(gainStepX, gainStepY)) * (j - (gainStepY - 1) * sideY)/sideY +...
39 (rGain(gainStepX+1, gainStepY+1) + rGain(gainStepX, gainStepY) - rGain(gainStepX+1, gainStepY)- rGain(gainStepX, gainStepY + 1)) *...
40 (i - (gainStepX - 1) * sideX)/sideX * (j - (gainStepY - 1) * sideY)/sideY + rGain(gainStepX, gainStepY);
41
42 gGainTmp = (gGain(gainStepX+1, gainStepY) - gGain(gainStepX, gainStepY)) * (i - (gainStepX - 1) * sideX)/sideX +...
43 (gGain(gainStepX, gainStepY+1) - gGain(gainStepX, gainStepY)) * (j - (gainStepY - 1) * sideY)/sideY +...
44 (gGain(gainStepX+1, gainStepY+1) + gGain(gainStepX, gainStepY) - gGain(gainStepX+1, gainStepY)- gGain(gainStepX, gainStepY + 1)) *...
45 (i - (gainStepX - 1) * sideX)/sideX * (j - (gainStepY - 1) * sideY)/sideY + gGain(gainStepX, gainStepY);
46
47 bGainTmp = (bGain(gainStepX+1, gainStepY) - bGain(gainStepX, gainStepY)) * (i - (gainStepX - 1) * sideX)/sideX +...
48 (bGain(gainStepX, gainStepY+1) - bGain(gainStepX, gainStepY)) * (j - (gainStepY - 1) * sideY)/sideY +...
49 (bGain(gainStepX+1, gainStepY+1) + bGain(gainStepX, gainStepY) - bGain(gainStepX+1, gainStepY)- rGain(gainStepX, gainStepY + 1)) *...
50 (i - (gainStepX - 1) * sideX)/sideX * (j - (gainStepY - 1) * sideY)/sideY + bGain(gainStepX, gainStepY);
51
52 disImg(i,j,1) = double(image_r(i, j)) * rGainTmp;
53 disImg(i,j,2) = double(image_g(i, j)) * gGainTmp;
54 disImg(i,j,3) = double(image_b(i, j)) * bGainTmp;
55 end
56end
57
58figure();
59subplot(121);imshow(image);title('org image');
60subplot(122);imshow(uint8(disImg));title('corrected image');