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

圖像處理中使用的導數

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

前言

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

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

　　設有定義域和取值都在實數域中的函數 y=f(x)y=f(x). 若 f(x)f(x) f(x); 在點 x0x0的某個鄰域內有定義，則當自變量 xx 在x0x0 處取得增量 ΔxΔx（點 x0+Δxx0+Δx 仍在該鄰域內）時，相應地 yy 取得增量Δy=f(x0+Δx)−f(x0)Δy=f(x0+Δx)−f(x0)；如果 ΔyΔy與 ΔxΔx 之比當 Δx→0Δx→0 時的極限存在，則稱函數 y=f(x)y=f(x) 在點 x0x0; 處可導，並稱這個極限為函數y=f(x)y=f(x) 在點 x0x0 處的導數，記為 f′(x0)f′(x0)，即：

 

f′(x0)=limΔx→0ΔyΔx=limΔx→0f(x0+Δx)−f(x0)Δxf′(x0)=limΔx→0ΔyΔx=limΔx→0f(x0+Δx)−f(x0)Δx

 

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

離散情況下的差分計算

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

前向差分求導　

df=x(n)−x(n−h)hdf=x(n)−x(n−h)h

后向差分求導

df=x(n+h)−x(n)hdf=x(n+h)−x(n)h

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

f(x+Δx)−f(x)Δx=−x−μσ2f(x)f(x+Δx)−f(x)Δx=−x−μσ2f(x)

　化簡可得：

f(x+Δx)=(1−x−μσ2Δx)f(x)f(x+Δx)=(1−x−μσ2Δx)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() 

執行結果：

圖像導數實現

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

 

gx=−1−1−1000111gx=−101−101−101

 

 

gy=−101−101−101gy=−1−1−1000111

 

下面我們用最后得出梯度的幅值為G(x,y)=(g2x+g2y)−−−−−−−−√G(x,y)=(gx2+gy2)方向為: θ=arctangygxθ=arctan⁡gygx現在我們用程序來實現這個過程。

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

 

L(x,y)=∂f(x)∂x(2)+∂f(y)∂y(2)L(x,y)=∂f(x)∂x(2)+∂f(y)∂y(2)

 

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

 

gxx=gx(x+1)−gx(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)gxx=gx(x+1)−gx(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)

 

最后同理可得：

 

gyy=f(x,y+1)−2∗f(x,y)+f(x,y−1)gyy=f(x,y+1)−2∗f(x,y)+f(x,y−1)

 

最后可得：

 

L(x,y)=∂f(x)∂x(2)+∂f(y)∂y(2)=f(x+1,y)−4∗f(x,y)+f(x−1,y)+f(x,y+1)+f(x,y−1)L(x,y)=∂f(x)∂x(2)+∂f(y)∂y(2)=f(x+1,y)−4∗f(x,y)+f(x−1,y)+f(x,y+1)+f(x,y−1)

 

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

 

0−10−14−10−1−00−10−14−10−1−0

 

最后代碼實現為：

""" @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矩陣，大體上按以上流程，基本都能實現計算出需要的模板。