isp算法学习-UVNR（UV color space Noise Reduce）
500 Words|Read in about 3 Min|本文总阅读量次
本文是isp算法学习的第十四篇章，主要是讲述UVNR（UV color space Noise Reduce）。
1简介
UVNR也就是经过CSC只有在YUV域对UV两个色域进行降噪，在有些方案里也叫CNR（chroma noise reduction）。主要就是在YUV域针对彩燥进行特殊处理的一系列算法。
关于噪声产生的原因在前面关于降噪的文章和视频中已经做了描述和讲解，这里就不做更多的解释，关于彩噪产生的原因，可以通过一个简单的方式来理解，就是如果一个灰色的画面，原本R，G，B三个通道的值应该是一样的（理论白平衡之后），但是由于前面提到过的噪声因素会导致三个通道的pixel value发生变化，且三个通道变化是随机的，那么三个通道各自偏差就不确定，那么经过插值处理后原本应该一样的三个通道就会出现偏差，也就出现了各种颜色，也就导致的彩燥的出现，然后经过pipeline各个阶段的gain的加持就会更显著。
Assuming a sensor with a color filter array (CFA) in a typical color imaging processing pipeline, the sensor gains and possibly other digital gains must be adjusted as a result of color correction, white balance and other steps which may account for unequal sensor sensitivity on the different channels. According to the propagation of error formula, signal variance increases with the square root of the applied gain, thus regardless of the analog or digital domain where such a gain adjustment is performed, the net effect is a different signal standard deviation on the various channels. After color interpolation, when sensor channels are correlated, the described mechanism gives rise to what it is referred to as chrominance noise. This is further amplified when a non-spatially adaptive color correction is applied, as noted by others in [1-4], to become large, perceptually periodic groupings of 15 to 25 pixels across. For completion and an in depth, low level discussion on the various physical sources of noise in an imaging system, the reader is referred to Janesik [5].
假设在一个典型的彩色成像处理管道中有一个带有彩色滤波器阵列(CFA)的传感器，由于色彩校正、白平衡和其他步骤，传感器增益和可能的其他数字增益必须进行调整，这可能导致不同通道上传感器灵敏度的不相等。根据误差传播公式，信号方差随施加的增益的平方根而增加，因此无论在模拟域还是数字域进行增益调整，其净效应是在各个通道上产生不同的信号标准差。**在颜色插值之后，当传感器通道相关时，所描述的机制会产生所谓的色度噪声。**当应用非空间自适应色彩校正时，这将进一步放大，正如其他人在[1-4]中所指出的那样，成为15到25像素的大型、感知周期性分组。为了完成对成像系统中各种物理噪声源的深入、低级讨论，请参考Janesik[5]。
2算法讲解
降噪算是图像处理一个很大的方向，所以这一块的算法很多，没法一一介绍，针对彩噪其实也就是将基本的降噪或者灰度图像降噪的方法进行调整后适配这个应用。这里就直接通过几个专利介绍几个相关算法。
2.1专利CHROMA NOISE REDUCTION FOR CAMERAS
这是苹果公司2009的一份专利，可以点击这里
算法主要思路如下，就是再YUV域对Cr和Cb通道分别进行处理。
比如针对Cb，在一个滤波窗口中，中心点的像素值就通过周围点的值的加权平均值来计算得到。权重的计算也很简单，就是判断该点和窗口中心的Cr,Cb插值的和是否大于阈值，如果大于阈值就认为是边缘，那么就不参与平均，权重设置为0，如果小于阈值就参与计算，权重就设置为1。
C'_{r_0}=\frac{\sum_{n=-N_1}^{N_2}(C_{r_n}T_n)}{\sum_{n=-N_1}^{N_2}(T_n)}\
$$
2.2专利IMAGE CHROMA NOISE REDUCTION
该专利是ST公司2012年的一份专利，，可以点击这里
首先计算出窗口中各个通道的分布范围：
$$
D_{Cr}=max(Cr)-min(Cr)\\
D_{Cb}=max(Cb)-min(Cb)\\
D_Y=max(Y)-min(Y)
$$
根据分布范围求出一个滤波强度系数：
$$
C_f=
\begin{cases}
D_Y,&&&if \ D_Y=min(D_Y,D_{Cr},D_{Cb})\\
max(D_Y,D_{Cr},D_{Cb})&&&otherwise
\end{cases}
$$
针对Cr和Cb**进行加权平均降噪**：
$$
denoisedCr=\sum\limits_{k = 1}^{M×M}\frac{WCr_k}{sumweightCr}·Cr[k]\\
denoisedCb=\sum\limits_{k = 1}^{M×M}\frac{WCb_k}{sumweightCb}·Cb[k]\\
$$
权重计算方式是通过亮度知道色度，然后通过高斯分布的方式计算，差别越大越不相似，权重就越小。
$$
WCr_k=e^{-\frac{1}{2}(\frac{|Y(k)-Y(c)|}{sigmaY})}·
e^{-\frac{1}{2}(\frac{|Cr(k)-Cr(c)|}{sigmaCr})}\\
sumweightCr=\sum\limits_{k = 1}^{M×M}WCr_k\\
WCb_k=e^{-\frac{1}{2}(\frac{|Y(k)-Y(c)|}{sigmaY})}·
e^{-\frac{1}{2}(\frac{|Cb(k)-Cb(c)|}{sigmaCb})}\\
sumweightCb=\sum\limits_{k = 1}^{M×M}WCb_k\\
$$
最后通过将原始图像和去噪后的图像进行blending得到最终的图像，系数由局部动态范围决定。
$$
Cr_k=originalCr_k+f(Cf)·(denoisedCr-originalCr_k)\\
Cb_k=originalCb_k+f(Cf)·(denoisedCb-originalCb_k)\\
f(x)=e^{\frac{1}{2}(\frac{x}{sigma})^2},\ x∈[0,max_value]
$$
2.3专利REMOVING CHROMA NOISE FROM DIGITAL IMAGES BY USIING VARIABLE SHAPE PIXEL NEIGHBORHOOD REGIONS
该专利是柯达公司1999年的一份专利
对应的算法
其实是一个均值滤波的变式，这个算法是作用到L*a*b*的颜色空间的，跟YUV还是有所区别的，但是万变不离其宗，分别对水平、竖直，两个斜45度进行运算。
$$
h=\frac{1}{13}
\begin{Bmatrix}
1&2&0&-2&-1\\
1&2&0&-2&-1\\
1&2&0&-2&-1\\
1&2&0&-2&-1\\
1&2&0&-2&-1
\end{Bmatrix}\\
v=\frac{1}{13}
\begin{Bmatrix}
-1&-1&-1&-1&-1\\
1&-2&-2&-2&-1\\
0&0&0&0&0\\
1&2&2&2&1\\
1&1&1&1&1
\end{Bmatrix}\\
s=\frac{1}{13}
\begin{Bmatrix}
-1&-1&-1&-1&0\\
1&-2&-2&0&1\\
-1&-2&0&2&1\\
-1&0&2&2&1\\
0&1&1&1&1
\end{Bmatrix}\\
b=\frac{1}{13}
\begin{Bmatrix}
0&-1&-1&-1&-1\\
1&0&-2&-2&1\\
1&2&0&-2&-1\\
1&2&2&0&-1\\
1&1&1&1&0
\end{Bmatrix}\\
$$
对应的公式为
$$
g(x)=|hf(x)|+|vf(x)|+|sf(x)|+|bf(x)|
$$
其中**指的是是二维卷积运算的绝对值
-
通过5*5的边缘滤波处理，区分通道
-
把各个通道的滤波结果相加求得edge map
-
得到整个图的Edge map后，开始遍历整个edge map
-
针对如图的具体的某个点A的处理方式
1）以A点为中心，沿着N，NE，SE，S，SW，W，NM八个方向进行检索
2）以A点在edge map中的value作为reference value
3）以检索点在edge map中的value与A的插值△作为判别依据
4）如果△<Th，那么就是平摊去，加入到待处理区域
如果△>Th，那么可能是边缘，就修正该方向的检索
5）当检索完成后就用该ROI中的值作为降噪后的值
3算法
1clc; close all; clear;
2
3%% ------- global value -----------
4TH = struct('ThL1', 0, 'ThL2', 0, 'ThL3', 0, 'ThL4', 0, 'ThL5', 0, 'ThL6', 0, 'ThL7', 0,
5 'ThH1', 0, 'ThH2', 0, 'ThH3', 0, 'ThH4', 0, 'ThH5', 0, 'ThH6', 0, 'ThH7', 0);
6KSharp = 128;
7kernelHSize = 3;
8kernelVSize = 3;
9GaussSigma = 0.5;
10
11%% ------- Data preparation -------
12imgRGB = imread('./images/kodim19.png');
13[h, w, c] = size(imgRGB);
14imgHSV = rgb2hsv(imgRGB);
15padHSize = floor(kernelHSize / 2);
16padVSize = floor(kernelVSize / 2);
17imgHSVPad = padarray(imgHSV, [padHSize, padVSize], 'replicate');
18dstHSV = zeros(h, w, c);
19guassKernel = fspecial('gaussian',[kernelVSize kernelHSize], GaussSigma);
20
21Eedge = zeros(h, w);
22%% ----------- main loop ----------
23for y = 1+padHSize: w+padHSize
24 for x = 1+padVSize: h+padVSize
25 calWin = imgHSVPad(x-padVSize: x+padVSize, y-padHSize: y+padHSize, :);
26 lp = calLp(calWin, TH);
27 hp = calHp(calWin, TH);
28 dstHSV(x-padHSize, y-padVSize, 1) = sum(sum(lp .* calWin(:, :, 1))) / sum(sum(lp));
29 dstHSV(x-padHSize, y-padVSize, 2) = sum(sum(lp .* calWin(:, :, 2))) / sum(sum(lp));
30 Eedge(x-padHSize, y-padVSize) = sum(sum(hp .* calWin(:, :, 2))) * KSharp/ sum(sum(hp));
31 end
32end
