Guided Filter

Posted on Edited on

输入图p，导向图I， 输出图q，\(\Rightarrow q_i = \sum_jW_{ij}(I)p_j\)，其中i，j为pixel index，\(W_{ij}\)是导向图I的一个函数，与输入图p独立。

假设输出图q在一个窗口\(w_i\)中是导向图I的相信变换，即：\(q_i = a_kI_i + b_k, \forall i \in w_k\)

这种模型保证了仅有在导向图I有边缘的时候，输出图也有边缘，\(\nabla q = a\nabla I\)，原理可以按照https://sci-hubtw.hkvisa.net/10.1109/acv.2002.1182150

为了确定\(a_k,b_k\)系数，有如下的优化函数。

\(E(a_k, b_k) = \sum_{i \in \omega_k}((a_kI_i + b_k - p_i)^2 + \epsilon a_k^2)\)

最终得到\(a_k,b_k\)如下：

\(a_k = \frac{\frac{1}{\lvert\epsilon\rvert}\sum{i\in\omega_k(I_ip_i - \mu_k\overline{p}_k)}}{\sigma_k^2 + \omega}\)

\(b_k = \overline{p_k} - a_k\mu_k\)

以偏导的方式求解\(a_k，b_k\)

对于输出图q，可以认为是输入图p减去其噪声，\(q_i = p_i - n_i\)，\(n_i\)表示噪声

对于输出图q与导向图I，有如下公式：\(\nabla q_i = a \nabla I_i\)， \(q_i = aI_i + b\)

最小化输出图q与输入图p的差异，求解系数：

\(\frac{\partial E}{\partial a_k} = 2\sum_{i}^{N}((a_kI_i + b_k - p_i)I_i + \epsilon a_k)\)

\(\frac{\partial E}{\partial b_k} = -2\sum_{i}^{N}(p_i - a_kI_i - b_k)\)

求解\(b_k\)

\[\begin{aligned} &\ \frac{\partial E}{\partial b_k} = -2\sum_{i}^{N}(p_i - a_kI_i - b_k) = 0 \\ &\ \Rightarrow Nb_k = \sum_{i}^n(p_i - a_kI_i) \\ &\ \Rightarrow b_k = \frac{1}{N} \sum_{i}^n(p_i - a_kI_i) \\ &\ \Rightarrow b_k = \frac{1}{N} \sum_{i}^np_i - a_k\frac{1}{N} \sum_{i}^nI_i \\ &\ \Rightarrow b_k = p_k - a_k u_k \\ \end{aligned}\]

\(p_k\)和\(u_k\)是窗口\(w_i\)内，输入图\(p_i\)和导向图\(I_i\)的均值。

求解\(a_k\)

\[\begin{aligned} &\ \frac{\partial E}{\partial a_k} = 2\sum_{i}^{N}((a_kI_i + b_k - p_i)I_i + \epsilon a_k) = 0 \\ &\ \Rightarrow -\sum_{i}^Np_iI_i + \sum_{i}^Nb_kI_i + \sum_{i}^Na_kI_i^2 + \sum_{i}^N \epsilon a_k = 0 \\ &\ \Rightarrow \sum_{i}^Np_iI_i - \sum_{i}^Nb_kI_i = \sum_{i}^Na_kI_i^2 + \sum_{i}^2\epsilon a_k \\ &\ \Rightarrow \sum_{i}^Np_iI_i - \sum_{i}^N(p_k-a_ku_k)I_i = a_k\sum_{i}^N(I_i^2 + \epsilon) \\ &\ \Rightarrow \sum_{i}^Np_iI_i - \sum_{i}^Np_kI_i + \sum_{i}^Na_ku_kI_i = a_k\sum_{i}^N(I_i^2 + \epsilon) \\ &\ \Rightarrow \sum_{i}^Np_iI_i - p_k\sum_{i}^NI_i = a_k\sum_{i}^N(I_i^2 + \epsilon - u_kI_i) \\ &\ \Rightarrow a_k = \frac{\sum_{i}^Np_iI_i - p_k\sum_{i}^NI_i}{\sum_{i}^N(I_i^2 + \epsilon - u_kI_i)} \\ &\ \Rightarrow a_k = \frac{\sum_{i}^Np_iI_i - \frac{1}{N}\sum_{i}^Np_i\sum_{i}^NI_i}{\sum_{i}^NI_i^2 - \frac{1}{N}\sum_{i}^NI_i\sum_{i}^NI_i + \epsilon} \\ &\ \Rightarrow a_k = \frac{\frac{1}{N}\sum_{i}^Np_iI_i - \frac{1}{N}\sum_{i}^Np_i\frac{1}{N}\sum_{i}^NI_i}{\frac{1}{N}\sum_{i}^NI_i^2 - \frac{1}{N}\sum_{i}^NI_i\frac{1}{N}\sum_{i}^NI_i + \epsilon} \\ \end{aligned}\]

根据协方差公式

\(var(X) = \frac{\sum_{i}^{N}(X_i-\overline X)(X_i - \overline X)}{n-1}\)

\(cov(X, Y) = E[(X - E[X])(Y - E[Y])]\)

\(cov(X, Y) = E[XY] - E[X][Y]\)

\(a_k = \frac{cov(p_i, I_i)}{\sigma_k^2 + \epsilon}\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
cv::Mat guidedFilter(cv::Mat I, cv::Mat p, int r, double eps)
{
  /*
  % GUIDEDFILTER   O(1) time implementation of guided filter.
  %
  % - guidance image: I (should be a gray-scale/single channel image)
  % - filtering input image: p (should be a gray-scale/single channel image)
  % - local window radius: r
  % - regularization parameter: eps
  */
  cv::Mat _I;
  I.convertTo(_I, CV_64FC1);
  I = _I;
  cv::Mat _p;
  p.convertTo(_p, CV_64FC1);
  p = _p;
  //[hei, wid] = size(I);
  int hei = I.rows;
  int wid = I.cols;
  //N = boxfilter(ones(hei, wid), r); % the size of each local patch; N=(2r+1)^2 except for boundary pixels.
  cv::Mat N;
  cv::boxFilter(cv::Mat::ones(hei, wid,I.type()), N, CV_64FC1, cv::Size(r, r));
  //mean_I = boxfilter(I, r) ./ N;
  cv::Mat mean_I;
  cv::boxFilter(I, mean_I, CV_64FC1, cv::Size(r,r));
  //mean_p = boxfilter(p, r) ./ N;
  cv::Mat mean_p;
  cv::boxFilter(p, mean_p, CV_64FC1, cv::Size(r,r));
  //mean_Ip = boxfilter(I.*p, r) ./ N;
  cv::Mat mean_Ip;
  cv::boxFilter(I.mul(p), mean_Ip, CV_64FC1,
  cv::Size(r, r));
  //cov_Ip = mean_Ip - mean_I .* mean_p; % this is the covariance of (I, p) in each local patch.
  cv::Mat cov_Ip = mean_Ip - mean_I.mul(mean_p);
  //mean_II = boxfilter(I.*I, r) ./ N;
  cv::Mat mean_II;
  cv::boxFilter(I.mul(I), mean_II, CV_64FC1,
  cv::Size(r, r));
  //var_I = mean_II - mean_I .* mean_I;
  cv::Mat var_I = mean_II - mean_I.mul(mean_I);
  //a = cov_Ip ./ (var_I + eps); % Eqn. (5) in the paper;
  cv::Mat a = cov_Ip / (var_I + eps);
  //b = mean_p - a .* mean_I; % Eqn. (6) in the paper;
  cv::Mat b = mean_p - a.mul(mean_I);
  //mean_a = boxfilter(a, r) ./ N;
  cv::Mat mean_a;
  cv::boxFilter(a, mean_a, CV_64FC1, cv::Size(r,r));
  mean_a = mean_a / N;
  //mean_b = boxfilter(b, r) ./ N;
  cv::Mat mean_b;
  cv::boxFilter(b, mean_b, CV_64FC1, cv::Size(r,r));
  mean_b = mean_b / N;
  //q = mean_a .* I + mean_b; % Eqn. (8) in the paper;
  cv::Mat q = mean_a.mul(I) + mean_b;
  return q;
}