圖像處理中導數和模板的求法

圖像處理中使用的導數

copyright 版權所有，嚴禁抄襲，轉載需獲得本人授權,郵箱：zhaogoodwell@gmail

前言

　　工欲善其事必先利其器，在圖像處理中最常用的數學基礎有導數、卷積。今天我們主要討論下數字圖像處理中的導數，從從連續函數的導數概念出發，再到離散情況下的導數，最后使用代碼來實現。所有只講理論，不給處實例代碼的行為都是耍流氓！！！

連續函數導數的一般性定義

　　設有定義域和取值都在實數域中的函數 \(y=f(x)\). 若 \(f(x)\) f(x); 在點 \(x_{0}\)的某個鄰域內有定義，則當自變量 \(x\) 在$ x_{0}$ 處取得增量 $\Delta x $（點 $ x_{0}+\Delta x$ 仍在該鄰域內）時，相應地 \(y\) 取得增量$ \Delta y=f(x_{0}+\Delta x)-f(x_{0})$；如果 \(\Delta y\)與 \(\Delta x\) 之比當 \(\Delta x \to 0\) 時的極限存在，則稱函數 \(y=f(x)\) 在點 \(x_{0}\); 處可導，並稱這個極限為函數\(y=f(x)\) 在點 \(x_{0}\) 處的導數，記為 \(f'(x_{0})\)，即：

\[f'(x_{0})=\lim _{\Delta x\to 0} {\frac {\Delta y}{\Delta x}} = \lim _{\Delta x\to 0}{\frac {f(x_{0}+\Delta x)-f(x_{0})}{\Delta x}} \]

這是導數的定義，需要用到極限，顯然這兒公式沒法在離散情況下套用，在離散情況下我們怎么來計算導數呢？差分。

離散情況下的差分計算

　　在離散情況下我們利用差分來代替微分，差分分為兩種，前向差分和后向差分。我們假設有一個數列\(x(n)\),\(n,h\in N^{+}\),我們有：

前向差分求導　$$df = \frac{x(n)-x(n-h)}{h}$$

后向差分求導 $$df = \frac{x(n+h) - x(n)}{h}$$

差分求導和導數形式上是一樣的，區別在於沒有使用極限，所以利用差分來求導是一種不精確的方法，我們把$h$稱作步長，顯然步長越小越好，就越接近導數的真實值，精度會越高。差分求導在數值分析里面有廣泛的使用。我們就拿一個微分方程仿真為例，來說明如何使用離散導數解決問題。我們知道$f(x) = \exp{\left(-\frac{(x-\mu)^2}{2\sigma ^2} \right)}$的導數是$f'(x) = -\frac{x-\mu}{\sigma ^2} f(x)$ 我們可以使用差分來仿真這個微分方程。對這個方程改寫如下：

$$\frac{f(x+\Delta x) -f(x)}{\Delta x} = -\frac{x-\mu}{\sigma ^2} f(x)$$

　化簡可得：

$$f(x+\Delta x) = \left(1 - \frac{x-\mu}{\sigma ^2} \Delta x \right) f(x)$$

　 我們可以設置一個很小的步長每次對進行迭代，可以得出函數的數值解。

我們使用如下Python 代碼對它進行仿真。

# -*- coding: utf-8 -*-
"""
@author: zhao
"""


import numpy as np
import matplotlib.pyplot as plt

def g(x):
    return np.exp(-((x-mu)**2)/(2*sigma**2))

sigma = 1.6
mu = 2
x = np.linspace(-6,10,200)


plt.plot(x,g(x))
plt.title('g(x)')
plt.show()



it = [200,2000,20000]
for step in it:
    x = np.linspace(-6,10,step)
    f = np.zeros(x.shape)
    delta = x[1] - x[0]
    f[0] = g(-6)
    for i in np.arange(1,step):
        f[i] = (1 - delta * (x[i-1]-mu)/(sigma**2)) * f[i-1]
    plt.plot(x,f)
    plt.title("step="+str(step))
    plt.show()
    

執行結果：

圖像導數實現

我們對圖像很多操作都是用模板來實現的，比如圖像的梯度，濾波，邊沿提取等技術。我們所說的圖像處理一般是指數字圖像，是對模擬信號的采樣，對圖像進行求導的操作只能通過差分等方式來實現。對待一幅圖像我們定義它的\(x\)方向上的導數為\(g_{x} = f(x+1) - f(x)\)，但是這樣沒有一個中心點我們操作起來不方便，所以我們就把這個模板進行擴展，所以我們采用如下模板來計算圖像的\(x\)和\(y\)方向的梯度:

\[g_{x} = \begin{array}{|c|c|c|} \hline -1& 0&1 \\ \hline -1& 0&1 \\ \hline -1& 0&1 \\ \hline \end{array} \]

\[g_{y} = \begin{array}{|c|c|c|} \hline -1& -1&-1 \\ \hline 0& 0&0 \\ \hline 1& 1&1 \\ \hline \end{array} \]

下面我們用最后得出梯度的幅值為\(G(x,y) = \sqrt{ \left(g_{x}^2 +g_{y}^2\right)}\)方向為: $ \theta = \arctan{\frac{g_{y}}{g_{x}}}$現在我們用程序來實現這個過程。

拉普拉斯算子，在數學上的表達式為：

\[L(x,y) = \frac{\partial f(x)}{\partial x^(2)} + \frac{\partial f(y)}{\partial y^(2)} \]

這個是對圖像\(x\)和\(y\)方向兩次求導，然后相加。我門先看\(x\)方向的一階導數，\(g_{x} = f(x,y) - f(x-1,y)\)，再對以一階導數求導便是二階導數,最終結果為：

\[g_{xx} = g_{x}(x+1) - g_{x}(x) = f(x+1,y)-f(x,y) - (f(x,y) - f(x-1,y)) = f(x+1,y) - 2*f(x,y) +f(x-1,y) \]

最后同理可得：

\[g_{yy} = f(x,y+1) - 2*f(x,y) +f(x,y -1) \]

最后可得：

\[L(x,y) = \frac{\partial f(x)}{\partial x^(2)} + \frac{\partial f(y)}{\partial y^(2)} = f(x+1,y) - 4*f(x,y) +f(x-1,y) + f(x,y+1) +f(x,y -1) \]

用3x的模板可以表示為：

\[\begin{array}{|c|c|c|} \hline 0& -1&0 \\ \hline -1& 4&-1 \\ \hline 0& -1&-0 \\ \hline \end{array} \]

最后代碼實現為：

"""
@author: zhao
"""

import numpy as np 
import cv2
import matplotlib.pyplot as plt

img = cv2.imread('lena1.tiff')
img = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
img = np.array(img,dtype= np.float64)
g_x = np.array([[-1,0,1],[-1,0,1],[-1,0,1]])
g_y = np.array([[-1,-1,-1],[0,0,0],[1,1,1]])
laplace = np.array([[0,-1,0],[-1,4,-1],[0,-1,0]])

img_g_x =cv2.filter2D(img,-1,g_x)
img_g_y =cv2.filter2D(img,-1,g_y)
img_laplace = cv2.filter2D(img,-1,laplace)

img_graid = np.sqrt(img_g_x **2 + img_g_y **2)
img_angle = np.arctan(img_g_y/(img_g_x + 2**-1000))

plt.subplot(2,3,1),plt.imshow(img,cmap ='gray'),plt.title('source'),plt.xticks([]),plt.yticks([])
plt.subplot(2,3,2),plt.imshow(img_g_x,cmap ='gray'),plt.title('x gradient'),plt.xticks([]),plt.yticks([])
plt.subplot(2,3,3),plt.imshow(img_g_y,cmap ='gray'),plt.title('y gradient'),plt.xticks([]),plt.yticks([])
plt.subplot(2,3,4),plt.imshow(img_graid,cmap ='gray'),plt.title('gradient  amplitude'),plt.xticks([]),plt.yticks([])
plt.subplot(2,3,5),plt.imshow(img_angle,cmap ='gray'),plt.title('angle'),plt.xticks([]),plt.yticks([])
plt.subplot(2,3,6),plt.imshow(img_laplace,cmap ='gray'),plt.title('Laplace'),plt.xticks([]),plt.yticks([])

在圖像處理中里面有很多跟導數有關的模板，比如在SIFT代碼中需要hessian矩陣，大體上按以上流程，基本都能實現計算出需要的模板。