拉普拉斯算子

在介绍拉普拉斯算子概念之前我们先介绍，哈密尔顿算子（$\nabla$），梯度，散度等概念

哈密尔顿算子

所谓哈密尔顿算子即为某一物理量在笛卡尔坐标系下的偏导数的矢量和，其运算符号为：$\nabla$，定义如下：
$$\nabla =\frac{\delta }{\delta x}i\text{i}i+\frac{\delta }{\delta y}j\text{j}j+\frac{\delta }{\delta z}k\text{k}k$$

梯度

当哈密尔顿算子$\nabla$作用于一个标量，比如函数$f$时，既可以表示为该标量在该空间上的梯度
简而言之，梯度即为$\nabla$作用于某一标量的结果，其表示如下：

$$gradf=\nabla f=\frac{\delta f}{\delta x_i}=\frac{\delta f}{\delta x}i\text{i}i+\frac{\delta f}{\delta y}j\text{j}j+\frac{\delta f}{\delta z}k\text{k}k$$

从上式不难发现，一个标量的梯度是一个矢量，它表示$f$在空间中某一点沿着某一方向的变化量，如果想要得到$f$在某一特定的方向$v$（即不是笛卡尔标系的方向）的梯度，即方向导数。
方向导数可以根据梯度与该方向$v$的点乘得到，其中$dt$表示沿着$v\text{v}v$方向的步长：

$$\frac{d{f}_v}{dt}=v\text{v}v\cdot \nabla f=\Vert \nabla f\Vert \cos (\nabla f,v\text{v}v)$$

由上式可得，当$\cos (\nabla f,v\text{v}v)=1$时，即空间上的某一方向与梯度方向一致时，表示沿该方向上有最大梯度，这也正符合我们所熟知的在空间上任一点的梯度为该点在空间上变化量最大的方向

散度

根据矢量的运算规则，当$\nabla$与一个矢量相乘时，其结果是一个标量，它代表着该矢量的散度：

$$divV=\nabla \cdot V=\frac{\delta V}{\delta x_i}=\frac{\delta u}{\delta x}+\frac{\delta v}{\delta y}+\frac{\delta w}{\delta z}$$

其中$V=u+v+w$表示一个矢量

关于散度大小定义：

$$div=\{\begin{array}{c}>0,向外发散\\ <0,向内聚合\\ =0,该点无源\end{array}$$

拉普拉斯算子

我们已经只带散度是由$\nabla$与一个矢量相乘时得到的，因此梯度作为一个矢量我们可以对其继续求散度，从而引入拉普拉斯算子${\nabla }^2$

$$\nabla \cdot \nabla f={\nabla }^{2}f=\frac{{\delta }^{2}f}{\delta x^2}+\frac{{\delta }^{2}f}{\delta y^2}+\frac{{\delta }^{2}f}{\delta z^2}$$

拉普拉斯算子从某种程度上看是对某一函数求二阶导数，我们在学数学分析的时候知道一阶导数表示函数某一点的变化量，而二阶导数则表示该点的变化速率，即凹凸性，我们将其放在图像上观察

右图图为左图的向量场，即由梯度组成的向量场，我们不难发现，向量场箭头汇聚的地方正好对应着三维图中凸点，向外发散的地方对应着凹点。

图像处理——拉普拉斯算子

假设一张图片的的某一列值看做函数$f(x)$，我们对其进行泰勒展开可以得到：

$$f(x+\Delta x)=f(x)+\Delta x\frac{\delta f(x)}{\delta x}+\frac{(\Delta x{)}^2}{2!}\frac{{\delta }^{2}f(x)}{\delta x^2}+\frac{(\Delta x{)}^3}{3!}\frac{{\delta }^{3}f(x)}{\delta x^3}+\cdots =\sum _{n=0}^{\infty }\frac{(\Delta x{)}^n}{n!}\frac{{\delta }^{n}f(x)}{\delta x^n}$$

上式中，$\Delta x$是$f$取样间隔。由于我们是对图像进行处理，所以这个间隔的单位是像素。所以我们对$x$前面取样有$\Delta x=1$，对$x$后面取样有$\Delta x=-1$，当$\Delta x=1$时，上式可化为：

$$f(x+1)=f(x)+\frac{\delta f(x)}{\delta x}+\frac{1}{2!}\frac{{\delta }^{2}f(x)}{\delta x^2}+\frac{1}{3!}\frac{{\delta }^{3}f(x)}{\delta x^3}+\cdots =\sum _{n=0}^{\infty }\frac{1}{n!}\frac{{\delta }^{n}f(x)}{\delta x^n}$$

同样当$\Delta x=-1$时有：
$$f(x-1)=f(x)-\frac{\delta f(x)}{\delta x}+\frac{1}{2!}\frac{{\delta }^{2}f(x)}{\delta x^2}-\frac{1}{3!}\frac{{\delta }^{3}f(x)}{\delta x^3}+\cdots =\sum _{n=0}^{\infty }\frac{(-1{)}^n}{n!}\frac{{\delta }^{n}f(x)}{\delta x^n}$$

对于一阶导数，我们只使用线性项，向前差分可以得到：
$$\frac{\delta f(x)}{\delta x}=f^{\prime }(x)=f(x+1)-f(x)$$

向后差分可以得到：
$$\frac{\delta f(x)}{\delta x}=f^{\prime }(x)=f(x)-f(x-1)$$

中心差分可的得到：
$$\frac{\delta f(x)}{\delta x}=f^{\prime }(x)=\frac{f(x+1)-f(x-1)}{2}$$

我们仅仅使用一阶导数展开，一般来说，用来表示导数的泰勒级数的项越多，近似就越精确。包含更多的项意味着在近似中使用更多的点，产生更小的误差。然而，对于相同数量的点来说，已经证明中心差分的误差较小。因此导数通常表示为中心差分。

基于中心差分的二阶导数$\frac{{\delta }^{2}f(x)}{\delta x^2}$可以由上面的$f(x+1)$泰勒展开与$f(x-1)$泰勒展开相加得到：

$$\frac{{\delta }^{2}f(x)}{\delta x^2}=(f^{\prime }(x){)}^{\prime }=f(x+1)-2f(x)+f(x-1)$$

上述只是一维上的二阶导数，而图像往往是二维坐标，我们可以写出双变量形式的一阶导与二阶导：

$$\frac{{\delta }^{2}f(x,y)}{\delta x^2}=f(x+1,y)-2f(x,y)+f(x-1,y)$$
$$\frac{{\delta }^{2}f(x,y)}{\delta y^2}=f(x,y+1)-2f(x,y)+f(x,y-1)$$

不难想到拉普拉斯算子应用到图像上的形式为：

$${\nabla }^{2}f(x,y)=\frac{{\delta }^{2}f}{\delta x^2}+\frac{{\delta }^{2}f}{\delta y^2}=f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)-4f(x,y)$$

表现在图像上如下：

上右图便是拉普拉斯卷积核，该卷积核只考虑到了某一像素点的上下左右的位置，我们可以将其扩展为考虑到周围像素点的卷积核

如上图的卷积核，使用在图像上我们可以检测到孤立点。
为了在对图像做拉普拉斯卷积后为了可以更加突出孤立点，我们可以设置一个阈值，大于该阈值的像素点赋值为1，小于阈值的赋值为0.

import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
imgfile = "./ImageProcessAndPractice/images/cameraman.jpg"
OriginalPic = np.array(Image.open(imgfile).convert('L'), dtype=np.int)
img = np.zeros((OriginalPic.shape[0]+2, OriginalPic.shape[1]+2), np.int)
# im = mpimg.imread(imgfile)
# print(im)
# plt.imshow(im)
########  制造遍历图像  ###################
for i in range(1, img.shape[0]-1):
    for j in range(1, img.shape[1]-1):
        img[i][j] = OriginalPic[i-1][j-1]

# 定义阈值函数
def threshold(x,v):
    if x > v:
        x = 255
    if x < v:
        x = 0
    return x

LaplacePic = np.zeros((OriginalPic.shape[0], OriginalPic.shape[1]), dtype=np.uint8)
kernel = [[1, 1, 1], [1, -8, 1], [1, 1, 1]]
for i in range(0, LaplacePic.shape[0]):
    for j in range(0, LaplacePic.shape[1]):
        LaplacePic[i][j] = threshold(abs(np.sum(np.multiply(kernel, img[i:i+3, j:j+3]))), 200)
fig = plt.figure(figsize=(12, 6))
plt.subplot(1,2,1)
plt.imshow(OriginalPic, cmap='gray')
plt.axis('off')
plt.subplot(1,2,2)
plt.imshow(LaplacePic, cmap='gray')
plt.axis('off')

其输出如下

可以清晰的看见，我们将图片经过拉普拉斯卷积后，可以显示出图片的轮廓