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

图像处理中使用的导数

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