【动手学计算机视觉】第三讲：图像预处理之图像分割
1
介绍
提到图像分割，主要包含两个方面：
非语义分割
语义分割
首先，介绍一下非语义分割。
非语义分割在图像分割中所占比重更高，目前算法也非常多，研究时间较长，而且算法也比较成熟，此类图像分割目前的算法主要有以下几种：
阈值分割
阈值分割是图像分割中应用最多的一类，该算法思想比较简单，给定输入图像一个特定阈值，如果这个阈值可以是灰度值，也可以是梯度值，如果大于这个阈值，则设定为前景像素值，如果小于这个阈值则设定为背景像素值。
阈值设置为100对图像进行分割：
区域分割
区域分割算法中比较有代表性的算法有两种：区域生长和区域分裂合并。
区域生长算法的核心思想是给定子区域一个种子像素，作为生长的起点，然后将种子像素周围邻域中与种子像素有相同或相似性质的像素(可以根据预先设定的规则，比如基于灰度差)合并到种子所在的区域中。
区域分裂合并基本上就是区域生长的逆过程，从整个图像出发，不断分裂得到各个子区域，然后再把前景区域合并，实现目标提取。
聚类
聚类是一个应用非常广泛的无监督学习算法，该算法在图像分割领域也有较多的应用。聚类的核心思想就是利用样本的相似性，把相似的像素点聚合成同一个子区域。
边缘分割
这是图像分割中较为成熟，而且较为常用的一类算法。边缘分割主要利用图像在边缘处灰度级会发生突变来对图像进行分割。常用的方法是利用差分求图像梯度，而在物体边缘处，梯度幅值会较大，所以可以利用梯度阈值进行分割，得到物体的边缘。对于阶跃状边缘，其位置对应一阶导数的极值点，对应二阶导数的过零点(零交叉点)。因此常用微分算子进行边缘检测。常用的一阶微分算子有Roberts算子、Prewitt算子和Sobel算子，二阶微分算子有Laplace算子和Kirsh算子等。由于边缘和噪声都是灰度不连续点，在频域均为高频分量，直接采用微分运算难以克服噪声的影响。因此用微分算子检测边缘前要对图像进行平滑滤波。LoG算子和Canny算子是具有平滑功能的二阶和一阶微分算子，边缘检测效果较好，因此Canny算子也是应用较多的一种边缘分割算法。
直方图
与前面提到的算法不同，直方图图像分割算法利用统计信息对图像进行分割。通过统计图像中的像素，得到图像的灰度直方图，然后在直方图的波峰和波谷是用于定位图像中的簇。
水平集
水平集方法最初由Osher和Sethian提出，目的是用于界面追踪。在90年代末期被广泛应用在各种图像领域。这一方法能够在隐式有效的应对曲线/曲面演化问题。基本思想是用一个符号函数表示演化中的轮廓（曲线或曲面），其中符号函数的零水平面对应于实际的轮廓。这样对应于轮廓运动方程，可以容易的导出隐式曲线/曲面的相似曲线流，当应用在零水平面上将会反映轮廓自身的演化。水平集方法具有许多优点：它是隐式的，参数自由的，提供了一种估计演化中的几何性质的直接方法，能够改变拓扑结构并且是本质的。
语义分割和非语义分割的共同之处都是要分割出图像中物体的边缘，但是二者也有本质的区别，用通俗的话介绍就是非语义分割只想提取物体的边缘，但是不关注目标的类别。而语义分割不仅要提取到边缘像素级别，还要知道这个目标是什么。因此，非语义分割是一种图像基础处理技术，而语义分割是一种机器视觉技术，难度也更大一些，目前比较成熟且应用广泛的语义分割算法有以下几种：
Grab cut
Mask R-CNN
U-Net
FCN
SegNet
由于篇幅有限，所以在这里就展开介绍语义分割，后期有时间会单独对某几个算法进行详细解析，本文主要介绍非语义分割算法，本文就以2015年UCLA提出的一种新型、高效的图像分割算法--相位拉伸变换为例，详细介绍一下，并从头到尾实现一遍。
2
相位拉伸变换
相位拉伸变换(Phase Stretch Transform, PST)，是UCLA JalaliLab于2015年提出的一种新型图像分割算法[Edge Detection in Digital Images Using Dispersive Phase Stretch Transform]，该算法主要有两个显著优点：
速度快
精度高
思想简单
实现容易
1Edge Detection in Digital Images Using Dispersive Phase Stretch Transform
2http://downloads.hindawi.com/journals/ijbi/2015/687819.pdfPST算法中，首先使用定位核对原始图像进行平滑，然后通过非线性频率相关（离散）相位操作，称为相位拉伸变换(PST)。 PST将2D相位函数应用于频域中的图像。施加到图像的相位量取决于频率;也就是说，较高的相位量被应用于图像的较高频率特征。由于图像边缘包含更高频率的特征，因此PST通过将更多相位应用于更高频率的特征来强调图像中的边缘信息。可以通过对PST输出相位图像进行阈值处理来提取图像边缘。在阈值处理之后，通过形态学操作进一步处理二值图像以找到图像边缘。思想主要包含三个步骤：
非线性相位离散化
阈值化处理
形态学运算
下面来详细介绍一下。
相位拉伸变换，核心就是一个公式，
其中为
3
编程实践
PST算法中最核心的就是公式(1)，编程实现可以一步一步来实现公式中的每个模块。
首先导入需要的模块，
1import os
2import numpy as np
3import mahotas as mh
4import matplotlib.pylab as plt
5import cv2
定义全局变量，
1L = 0.5
2S = 0.48
3W= 12.14
4Threshold_min = -1
5Threshold_max = 0.0019
6FLAG = 1
计算公式中的核心参数r和theta，
1def cart2pol(x, y):
2 theta = np.arctan2(y, x)
3 rho = np.hypot(x, y)
4 return theta, rho
生成变量p和q,
1x = np.linspace(-L, L, img.shape[0])
2y = np.linspace(-L, L, img.shape[1])
3X, Y = np.meshgrid(x, y)
4p, q = X.T, y.T
5theta, rho = cart2pol(p, q)
接下来对公式(1)从右至左依次实现，
对输入图像进行快速傅里叶变换,
1orig = np.fft.fft2(img)
实现平滑滤波核L,
1expo = np.fft.fftshift(np.exp(-np.power((np.divide(rho, math.sqrt((LPF ** 2) / np.log(2)))), 2)))
对图像进行平滑处理，
1orig_filtered = np.real(np.fft.ifft2((np.multiply(orig, expo))))
实现相位核，
1PST_Kernel_1 = np.multiply(np.dot(rho, W), np.arctan(np.dot(rho, W))) - 0.5 * np.log(1 + np.power(np.dot(rho, W), 2))
2PST_Kernel = PST_Kernel_1 / np.max(PST_Kernel_1) * S
将前面实现的部分与相位核做乘积，
1temp = np.multiply(np.fft.fftshift(np.exp(-1j * PST_Kernel)), np.fft.fft2(orig_filtered))
对图像进行逆快速傅里叶变换，
1temp = np.multiply(np.fft.fftshift(np.exp(-1j * PST_Kernel)), np.fft.fft2(Image_orig_filtered))
2orig_filtered_PST = np.fft.ifft2(temp)
进行角运算，得到变换图像的相位，
1PHI_features = np.angle(Image_orig_filtered_PST)
对图像进行阈值化处理，
1features = np.zeros((PHI_features.shape[0], PHI_features.shape[1]))
2features[PHI_features > Threshold_max] = 1
3features[PHI_features < Threshold_min] = 1
4features[I < (np.amax(I) / 20)] = 0
应用二进制形态学操作来清除转换后的图像,
1out = features
2out = mh.thin(out, 1)
3out = mh.bwperim(out, 4)
4out = mh.thin(out, 1)
5out = mh.erode(out, np.ones((1, 1)))
到这里就完成了相位拉伸变换的核心部分，
1def phase_stretch_transform(img, LPF, S, W, threshold_min, threshold_max, flag):
2 L = 0.5
3 x = np.linspace(-L, L, img.shape[0])
4 y = np.linspace(-L, L, img.shape[1])
5 [X1, Y1] = (np.meshgrid(x, y))
6 X = X1.T
7 Y = Y1.T
8 theta, rho = cart2pol(X, Y)
9 orig = ((np.fft.fft2(img)))
10 expo = np.fft.fftshift(np.exp(-np.power((np.divide(rho, math.sqrt((LPF ** 2) / np.log(2)))), 2)))
11 orig_filtered = np.real(np.fft.ifft2((np.multiply(orig, expo))))
12 PST_Kernel_1 = np.multiply(np.dot(rho, W), np.arctan(np.dot(rho, W))) - 0.5 * np.log(
13 1 + np.power(np.dot(rho, W), 2))
14 PST_Kernel = PST_Kernel_1 / np.max(PST_Kernel_1) * S
15 temp = np.multiply(np.fft.fftshift(np.exp(-1j * PST_Kernel)), np.fft.fft2(orig_filtered))
16 orig_filtered_PST = np.fft.ifft2(temp)
17 PHI_features = np.angle(orig_filtered_PST)
18 if flag == 0:
19 out = PHI_features
20 else:
21 features = np.zeros((PHI_features.shape[0], PHI_features.shape[1]))
22 features[PHI_features > threshold_max] = 1
23 features[PHI_features < threshold_min] = 1
24 features[img < (np.amax(img) / 20)] = 0
25
26 out = features
27 out = mh.thin(out, 1)
28 out = mh.bwperim(out, 4)
29 out = mh.thin(out, 1)
30 out = mh.erode(out, np.ones((1, 1)))
31 return out, PST_Kernel
下面完成调用部分的功能，
首先读取函数并把图像转化为灰度图，
1Image_orig = mh.imread("./cameraman.tif")
2if Image_orig.ndim == 3:
3 Image_orig_grey = mh.colors.rgb2grey(Image_orig)
4else:
5 Image_orig_grey = Image_orig
调用前面的函数，对图像进行相位拉伸变换，
1edge, kernel = phase_stretch_transform(Image_orig_grey, LPF, Phase_strength, Warp_strength, Threshold_min, Threshold_max, Morph_flag)
显示图像，
1Overlay = mh.overlay(Image_orig_grey, edge)
2edge = edge.astype(np.uint8)*255
3plt.imshow(Edge)
4plt.show()
主函数的完整内容为，
1def main():
2 Image_orig = mh.imread("./cameraman.tif")
3 if Image_orig.ndim == 3:
4 Image_orig_grey = mh.colors.rgb2grey(Image_orig)
5 else:
6 Image_orig_grey = Image_orig
7 edge, kernel = phase_stretch_transform(Image_orig_grey, LPF, S, W, Threshold_min,
8 Threshold_max, FLAG)
9 Overlay = mh.overlay(Image_orig_grey, Edge)
10 Edge = Edge.astype(np.uint8)*255
11 plt.imshow(Edge)
12 plt.show()
分割结果为：
来都来了，点个在看再走吧~~~