图像处理中导数和模板的求法

图像处理中使用的导数

前言

　　工欲善其事必先利其器，在图像处理中最常用的数学基础有导数、卷积。今天我们主要讨论下数字图像处理中的导数，从从连续函数的导数概念出发，再到离散情况下的导数，最后使用代码来实现。所有只讲理论，不给处实例代码的行为都是耍流氓！！！

连续函数导数的一般性定义

　　设有定义域和取值都在实数域中的函数 (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 o 0) 时的极限存在，则称函数 (y=f(x)) 在点 (x_{0}); 处可导，并称这个极限为函数(y=f(x)) 在点 (x_{0}) 处的导数，记为 (f'(x_{0}))，即：

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

这是导数的定义，需要用到极限，显然这儿公式没法在离散情况下套用，在离散情况下我们怎么来计算导数呢？差分。

离散情况下的差分计算

　　在离散情况下我们利用差分来代替微分，差分分为两种，前向差分和后向差分。我们假设有一个数列(x(n)),(n,hin 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}{2sigma ^2} ight)}$的导数是$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 ight) 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} = egin{array}{|c|c|c|} hline -1& 0&1 \ hline -1& 0&1 \ hline -1& 0&1 \ hline end{array} ]

[g_{y} = egin{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 ight)})方向为: $ heta = 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的模板可以表示为：

[egin{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矩阵，大体上按以上流程，基本都能实现计算出需要的模板。