導向濾波公式推導和擴展

何凱明大神的代表作之一
論文地址：Guided Image Filtering

導向濾波的一般表達方式

$q_{i}=\sum_{j} W_{i j}(I) p_{j}$qi​=j∑​Wij​(I)pj​
其中 $q$q表示輸出， $p$p表示輸入， $I$I表示導向圖。

先驗假設

假設在局部範圍內，輸出圖與導向圖的關係可以用一個線性模型表示：
$q_{i}=a_{k} I_{i}+b_{k}, \forall i \in \omega_{k}$qi​=ak​Ii​+bk​,∀i∈ωk​
另外輸出圖是由輸入圖減去噪聲（需要被濾掉的部分）得到
$q_{i}=p_{i}-n_{i}$qi​=pi​−ni​

能量函數和求解

我們需要做的就是最小化能量函數
$E\left(a_{k}, b_{k}\right)=\sum_{i \in \omega_{k}}\left(\left(a_{k} I_{i}+b_{k}-p_{i}\right)^{2}+\epsilon a_{k}^{2}\right)$E(ak​,bk​)=i∈ωk​∑​((ak​Ii​+bk​−pi​)2+ϵak2​)
其中 $\epsilon$ϵ爲正則項，防止係數 $a_k$ak​過大。
對 $a_k$ak​和 $b_k$bk​求導：
$\frac{\partial E}{\partial a_{k}}=2 \sum_{i}^{N}\left(\left(-p_{i}+a_{k} I_{i}+b_{k}\right) I_{i}+\epsilon a_{k}\right)$∂ak​∂E​=2i∑N​((−pi​+ak​Ii​+bk​)Ii​+ϵak​)
$\frac{\partial E}{\partial b_{k}}=-2 \sum_{i}^{N}\left(p_{i}-a_{k} I_{i}-b_{k}\right)$∂bk​∂E​=−2i∑N​(pi​−ak​Ii​−bk​)
先求取 $b_k$bk​，令偏導數爲零：
$\begin{aligned} &0= \sum_{i}^{N}\left(p_{i}-a_{k} I_{i}-b_{k}\right)\\ &\Rightarrow N b_{k}=\sum_{i}^{N}\left(p_{i}-a_{k} I_{i}\right)\\ &\Rightarrow b_{k}=\frac{1}{N} \sum_{i}^{N}\left(p_{i}-a_{k} I_{i}\right)\\ &\Rightarrow b_{k}=\frac{1}{N} \sum_{i}^{N} p_{i}-a_{k} \frac{1}{N} \sum_{i}^{N} I_{i} \\ &\Rightarrow b_{k}= \bar{p_{k}}-a_{k} \mu_{k} \end{aligned}$​0=i∑N​(pi​−ak​Ii​−bk​)⇒Nbk​=i∑N​(pi​−ak​Ii​)⇒bk​=N1​i∑N​(pi​−ak​Ii​)⇒bk​=N1​i∑N​pi​−ak​N1​i∑N​Ii​⇒bk​=pk​ˉ​−ak​μk​​
其中 $\bar{p_k}$pk​ˉ​爲輸入圖窗口內的平均值， $\mu_k$μk​爲導向圖在窗口內的平均值。
同理，再來求 $a_k$ak​
$\begin{aligned} &0= \sum_{i}^{N}\left(\left(-p_{i}+a_{k} I_{i}+b_{k}\right) I_{i}+\epsilon a_{k}\right)\\ &\Rightarrow 0=\sum-p_{i} I_{i}+\sum a_{k} I_{i}^{2}+\sum b_{k} I_{i}+\sum \varepsilon a_{k}\\ &\Rightarrow \sum p_{i} I_{i}-\sum b_{k} I_{i}=a_{k} \sum\left(I_{i}^{2}+\epsilon\right)\\ &\Rightarrow \sum p_{i} I_{i}-\sum\left(p_{k}-a_{k} \mu_{k}\right) I_{i}=a_{k} \sum\left(I_{i}^{2}+\epsilon\right)\\ &\Rightarrow a_{k}=\frac{\sum\left(p_{i} I_{i}-p_{k} I_{i}\right)}{\sum\left(I_{i}^{2}-\mu_{k} I_{i}+\epsilon\right)}\\ &\Rightarrow a_{k}=\frac{\sum\left(p_{i} I_{i}\right)-\frac{1}{N} \sum p_{i} \sum I_{i}}{\sum I_{i}^{2}-\frac{1}{N} \sum I_{i} \sum I_{i}+\epsilon}=\frac{\frac{1}{N} \sum\left(p_{i} I_{i}\right)-\frac{1}{N} \sum p_{i} \frac{1}{N} \sum I_{i}}{\frac{1}{N} \sum I_{i}^{2}-\frac{1}{N} \sum I_{i} \frac{1}{N} \sum I_{i}+\epsilon} \end{aligned}$​0=i∑N​((−pi​+ak​Ii​+bk​)Ii​+ϵak​)⇒0=∑−pi​Ii​+∑ak​Ii2​+∑bk​Ii​+∑εak​⇒∑pi​Ii​−∑bk​Ii​=ak​∑(Ii2​+ϵ)⇒∑pi​Ii​−∑(pk​−ak​μk​)Ii​=ak​∑(Ii2​+ϵ)⇒ak​=∑(Ii2​−μk​Ii​+ϵ)∑(pi​Ii​−pk​Ii​)​⇒ak​=∑Ii2​−N1​∑Ii​∑Ii​+ϵ∑(pi​Ii​)−N1​∑pi​∑Ii​​=N1​∑Ii2​−N1​∑Ii​N1​∑Ii​+ϵN1​∑(pi​Ii​)−N1​∑pi​N1​∑Ii​​​
根據方差和協方差公式
$\begin{aligned} &\operatorname{var}(X)=\frac{\sum_{i=1}^{n}\left(X_{i}-\bar{X}\right)^2}{n-1}\\ &\operatorname{Cov}(X, Y)=E[(X-E[X])(Y-E[Y])]\\ &\operatorname{Cov}(X, Y)=E[X Y]-E[X] E[Y] \end{aligned}$​var(X)=n−1∑i=1n​(Xi​−Xˉ)2​Cov(X,Y)=E[(X−E[X])(Y−E[Y])]Cov(X,Y)=E[XY]−E[X]E[Y]​
可以得到
$a_{k}=\frac{\operatorname{Cov}(p, I)}{\operatorname{Var}(I)+N \epsilon}=\frac{\operatorname{Cov}(p, I)}{\sigma_{k}^{2}+\epsilon}$ak​=Var(I)+NϵCov(p,I)​=σk2​+ϵCov(p,I)​
當導向圖與輸入圖相同時，導向濾波就變成了以輸入圖自身爲導向的保邊濾波。上面表達式就會變成
$\begin{aligned} a_{k} &=\frac{\sigma_{k}^{2}}{\sigma_{k}^{2}+\varepsilon} \\ b_{k} &=\left(1-a_{k}\right) \mu_{k} \end{aligned}$ak​bk​​=σk2​+εσk2​​=(1−ak​)μk​​

• 當遇到邊緣區域，方差較大，此時 $\sigma_{k}^2>>\epsilon$σk2​>>ϵ，有 $q=p$q=p
• 反之，當方差較小， $\sigma_{k}^2<<\epsilon$σk2​<<ϵ，輸出均值 $q=\mu_k$q=μk​

其實更仔細思考一下，這裏對「邊緣」和「平坦」區域的定義實際上是取決於 $\epsilon$ϵ的取值，當區域方差遠大於 $\epsilon$ϵ時，則「定義爲」這次濾波中的邊緣區域，在濾波時會被保留，反之則是平坦區域。

思考

圖像濾波與線性迴歸的關係

線性迴歸：有一些離散點 $x_i$xi​，需要擬合一條直線 $y_i$yi​使得 $\epsilon=\sum_{i}(y_i-(ax_i+b))^2$ϵ=∑i​(yi​−(axi​+b))2最小。
圖像去噪：有一些含有噪聲的像素 $p_i$pi​，需要擬合一個線性函數找到沒有噪聲的圖像 $q_i$qi​，使得 $\epsilon=\sum_{i}(p_i-q_i)^2$ϵ=∑i​(pi​−qi​)2最小。

擴展

簡化計算

可以看到按照濾波公式，對於輸出的每一個像素，都需要計算輸入圖和導向圖在窗口內的均值，以此得到 $a_k$ak​和 $b_k$bk​。但是這其中有很多重複計算，此時可以通過計算經過所有覆蓋當前像素的窗口的平均 $a_k$ak​和 $b_k$bk​來計算單個像素的輸出
$q_{i}=\bar{a}_{i} I_{i}+\bar{b}_{i}$qi​=aˉi​Ii​+bˉi​
經過這個改變之後可能無法滿足最初的 $\nabla q = a\nabla I$∇q=a∇I，但是仍然可以保證在強邊緣附近導向圖的信息是可以被保留下來的 $\nabla q \approx \bar{a} \nabla I$∇q≈aˉ∇I。

濾波覈計算

根據最開始的一般表達式，可以知道輸出 $q_i$qi​是線性依賴於輸入 $p_i$pi​，所以要求取濾波核只需要使用 $q_i$qi​對 $p_i$pi​求導就可以。將上面表達式中的 $b_k$bk​替換掉，可以得到
$q_{i}=\frac{1}{|\omega|} \sum_{k \in \omega_{i}}\left(a_{k}\left(I_{i}-\mu_{k}\right)+\bar{p}_{k}\right)$qi​=∣ω∣1​k∈ωi​∑​(ak​(Ii​−μk​)+pˉ​k​)
$\frac{\partial q_{i}}{\partial p_{j}}=\frac{1}{|\omega|} \sum_{k \in \omega_{i}}\left(\frac{\partial a_{k}}{\partial p_{j}}\left(I_{i}-\mu_{k}\right)+\frac{\partial \bar{p}_{k}}{\partial p_{j}}\right)$∂pj​∂qi​​=∣ω∣1​k∈ωi​∑​(∂pj​∂ak​​(Ii​−μk​)+∂pj​∂pˉ​k​​)
其中
$\frac{\partial \bar{p}_{k}}{\partial p_{j}}=\frac{1}{|\omega|} \delta_{j \in \omega_{k}}=\frac{1}{|\omega|} \delta_{k \in \omega_{j}}$∂pj​∂pˉ​k​​=∣ω∣1​δj∈ωk​​=∣ω∣1​δk∈ωj​​
是一個狄拉克函數，在窗口內部爲1，其餘地方爲0。
在根據 $a_k$ak​的表達式對 $p_j$pj​求導
$\begin{aligned} \frac{\partial a_{k}}{\partial p_{j}} &=\frac{1}{\sigma_{k}^{2}+\epsilon}\left(\frac{1}{|\omega|} \sum_{i \in \omega_{k}} \frac{\partial p_{i}}{\partial p_{j}} I_{i}-\frac{\partial \bar{p}_{k}}{\partial p_{j}} \mu_{k}\right) \\ &=\frac{1}{\sigma_{k}^{2}+\epsilon}\left(\frac{1}{|\omega|} I_{j}-\frac{1}{|\omega|} \mu_{k}\right) \delta_{k \in \omega_{j}} \end{aligned}$∂pj​∂ak​​​=σk2​+ϵ1​(∣ω∣1​i∈ωk​∑​∂pj​∂pi​​Ii​−∂pj​∂pˉ​k​​μk​)=σk2​+ϵ1​(∣ω∣1​Ij​−∣ω∣1​μk​)δk∈ωj​​​
代回原式就得到了濾波核的表達式
$\frac{\partial q_{i}}{\partial p_{j}}=\frac{1}{|\omega|^{2}} \sum_{k \in \omega_{i}, k \in \omega_{j}}\left(1+\frac{\left(I_{i}-\mu_{k}\right)\left(I_{j}-\mu_{k}\right)}{\sigma_{k}^{2}+\epsilon}\right)$∂pj​∂qi​​=∣ω∣21​k∈ωi​,k∈ωj​∑​(1+σk2​+ϵ(Ii​−μk​)(Ij​−μk​)​)
根據濾波核表達式，當像素 $i$i和 $j$j都在窗口內部，但處於一個邊緣的同側時（ ) \frac{\partial q_{i}}{\partial p_{j}}=\frac{1}{|\omega|^{2}} \sum_{k \in \omega_{i}, k \in \omega_{j}}\left(1+\frac{\left(I_{i}-\mu_{k}\right)\left(I_{j}-\mu_{k}\right)}{\sigma_{k}^{2}+\epsilon}\right) ∂pj​∂qi​​=∣ω∣21​k∈ωi​,k∈ωj​∑​(1+σk2​+ϵ(Ii​−μk​)(Ij​−μk​)​)
根據濾波核表達式，當像素 $i$i和 $j$j都在窗口內部，但處於一個邊緣的同側時（顏色相近），括號內後半部分的計算結果爲正，融合權重很高；反之如果處於邊緣的異側，括號內後半部分結果爲負，融合權重就會小很多。

高斯導向濾波

由於導向濾波算法採用的是均值濾波，在紋理不那麼強或是均勻紋理時，導向濾波會退化成兩個均值濾波的串聯，多次均值濾波的串聯雖然可以近似高斯濾波，但是僅兩次均值濾波還是會造成權重在x和y方向附近比其他方向稍大，如圖
這是可以在能量函數前面加一個高斯核，讓權重可以更加均勻分佈
$E(</$