Python语言OpenCV开发之使用OpenCV处理图像(下)

前言

本章内容继续上一节，因为内容比较多，所以分成了两节来讲述，都是和图像处理相关的，本节主要学习图像直方图，图像变换，模板匹配，分水岭算法等知识。

正文

1、 图像直方图
通过图像直方图可以对整幅图像的灰度值有一个整体的了解，直方图的x轴是灰度值(0到255),y轴是图片中具有相同灰度值的点的数目。直方图就是对图像的另一种解释。通过直方图可以对图像的对比度，亮度，灰度分布等有一个直观的认识。
1.1 直方图的统计
OpenCV统计直方图的函数是cv2.calcHist(images, channels, mask, histSize, ranges[, hist[, accumuate]])
参数images：原图像(图像格式为uint8或float32)。当传入函数时应该用中括号[]括起来。
参数channels：同样需要用中括号括起来，它会告诉函数要统计哪副图像的直方图。如果输入的图像是灰度图，它的值就是[0],如果是彩色图像的话，传入的参数可以是[0],[1],[2]分别对应着通道B,G,R。
参数mask：掩模图像。要统计整幅图像的直方图就把它设为None。但是要想统计图像某一部分的直方图的话，就需要制作一个掩模图像，并使用它。
参数histSize：BIN的数目。也应该用[]括起来。
参数ranges：像素值的范围，通常为[0,256]。
例如：

img = cv2.imread("xxx.jpg", 0)
hist = cv2.calcHist([img], [0], None, [256], [0,256])

当然也可以使用numpy中的方法做直方图的统计，函数是np.histogram()；

1.2 绘制直方图
绘制直方图有两种方法，简单点的就是使用Matplotlib中的绘制函数；例如：

img = cv2.imread('image/apple.jpg', 0)
plt.hist(img.ravel(), 256, [0,256])
plt.show()

还可以使用OpenCV中的绘制函数。例如：

img = cv2.imread('image/apple.jpg')
color = ('b', 'g', 'r')
for i, col in enumerate(color):
    histr = cv2.calcHist([img], [i], None, [256], [0,256])
    plt.plot(histr, color=col)
    plt.xlim([0,256])
plt.show()

1.3 使用掩模
要统计图像某个局部区域的直方图只需要构建一副掩模图像。将要统计的部分设置成白色，其余部分为黑色，就构成了一副掩模图像，传给函数即可；例如：

img = cv2.imread('image/apple.jpg', 0)
mask = np.zeros(img.shape[:2], np.uint8)
mask[40:180, 30:170] = 255
masked_img = cv2.bitwise_and(img, img, mask=mask)
hist_full = cv2.calcHist([img], [0], None, [256], [0,256])
hist_mask = cv2.calcHist([img], [0], mask, [256], [0,256])
plt.subplot(221),plt.imshow(img, 'gray')
plt.subplot(222),plt.imshow(mask, 'gray')
plt.subplot(223),plt.imshow(masked_img, 'gray')
plt.subplot(224),plt.plot(hist_full),plt.plot(hist_mask)
plt.xlim([0,256])
# plt.hist(img.ravel(), 256, [0,256])
plt.show()

1.4 直方图均衡化
如果一副图片整体很亮，那么所有的像素值应该都会很高。但是一副高质量的图像的像素值分布应该很广泛，so应该把它的直方图做一个横向拉伸，这就是直方图均衡化要做的事情。通常情况下这种操作会改善图像的对比度。就是用一个变换函数帮助我们把现有的直方图映射到一个广泛分布的直方图中，这就是直方图均衡化要做的事情。
即使输入的图片是一个比较暗的图片，在经过直方图均衡化之后也能得到相同的结果。因此直方图均衡化经常用来使所有的图片具有相同的亮度条件的参考工具。这在很多情况下都很有用，例如，面部识别，训练分类器前，训练集的所有图片都要先进行直方图均衡化从而使它们达到相同的亮度条件。
在OpenCV中使用函数：cv2.equalizeHist()可对图像进行均衡化。这个函数的输入图片仅仅是一副灰度图像，输出结果是直方图均衡化之后的图像。当然这是对一整副图像进行均衡化，还有一种方法是自适应的直方图均衡化，效果更好。这种情况下，整幅图像会被分成很多小块，这些小块被称为“tiles”（在 OpenCV 中 tiles 的大小默认是 8x8），然后再对每一个小块分别进行直方图均衡化（跟前面类似）。所以在每一个的区域中，直方图会集中在某一个小的区域中（除非有噪声干扰）。如果有噪声的话，噪声会被放大。为了避免这种情况的出现要使用对比度限制。对于每个小块来说，如果直方图中的 bin 超过对比度的上限的话，就把其中的像素点均匀分散到其他 bins 中，然后在进行直方图均衡化。最后，为了去除每一个小块之间“人造的”（由于算法造成）边界，再使用双线性差值，对小块进行缝合。实例：

ret1 = cv2.createCLAHE(clipLimit=2.0, titleGridSize=(8,8))
ret = ret1.apply(img)

1.5 2D直方图
前面绘制的是一维直方图，之所以称为一维，是因为只考虑了图像的一个特征：灰度值。2D直方图要考虑两个图像特征。对于彩色图像的直方图通常要考虑每个颜色的Hue和饱和度(Stauration)。
使用函数：cv2.calcHist()来计算直方图既简单又方便。如果要绘制颜色直方图的话，首先将图像的颜色空间从BGR转换到HSV（记住，计算一维直方图，要从BGR转换到HSV）。计算2D直方图，函数的参数要做如下：参数channels=[0, 1] 因为需要同时处理H和S两个通道；参数bins=[180,256] H通道为180，S通道为256；参数range=[0,180, 0,256]H的取值范围在0到180，S的取值范围在0到256；示例如下：

img = cv2.imread('image/apple.jpg')
hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)
hist = cv2.calcHist([hsv], [0,1], None, [180,256], [0,180,0,256])
plt.imshow(hist)

X轴显示S值，Y轴显示H值

1.6 直方图反向投影
可以用来做图像分割，或者在图像中找寻我们感兴趣的部分。简单来说，他会输出和输入图像同样大小的图像，其中每个像素值代表图像上对应点属于目标对象的概率。用更简单的话解释：输出图像中像素值越高（越白）的点就越可能代表我们要搜索的目标。
OpenCV提供的函数cv2.calcBackProject()可以用来做直方图反向投影。它的参数与函数cv2.calcHist的参数基本相同。其中的一个参数是我们要查找目标的直方图。同样再使用目标的直方图做反向投影之前我们应该先对其做归一化处理。返回的结果是一个概率图像，我们使用一个圆盘形卷积核对其做卷积操作，最后使用阈值进行二值化。

2. 图像变换
也可以说是傅里叶变换，经常用来分析不同滤波器的频率特性。可以使用2D离散傅里叶变换(DFT)分析图像的频域特性。实现DFT的一个快速算法称为快速傅里叶变换(FFT)。
对于一个正弦信号：x(t) = Asin(2πft),它的频率为f，如果把这个信号转到它的频域表示，我们会在频率f中看到一个峰值。如果我们的信号是由采样产生的离散信号组成的，我们会得到类似的频谱图，只不过前面是连续的，现在是离散的。你可以把图像想象成沿着两个方向采集的信号。所以对图像同时进行X方向和Y方向的傅里叶变换，就会得到这幅图像的频域表示（频谱图）。
关于傅里叶变换准备介绍两种方法：Numpy中的和OpenCV中的方法
2.1 Numpy中的傅里叶变换
Numpy中的FFT包可以实现快速傅里叶变换。函数np.fft.fft2()可以对信号进行频率转换，输出结果是一个复杂的数组。第一个参数是输入图像，要求是灰度格式。第二个参数是可选的，决定输出数组的大小。输出数组的大小和输入图像大小一样。如果输出结果比输入图像大，输入图像就需要进行FFT前补0.如果输出结果比输入图像小的话，输入图像就会被切割。
现在我们得到了结果，频率为0的部分（直流分量）在输出图像的左上角，如果想让它（直流分量）在输出图像的中心，还需要将结果沿两个方向平移N/2。函数np.fft.fftshift()可以实现，进行完频率变换后，就可以构建振幅谱了；实例如下：

import cv2 
import numpy as np 
from matplotlib import pyplot as plt
img = cv2.imread('image/rose.jpg', 0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
# 这里构建振幅图
magnitude_spectrum = 20*np.log(np.abs(fshift))
plt.subplot(121),plt.imshow(img, cmap='gray')
plt.title("Input Image"),plt.xticks([]),plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap='gray')
plt.title("Magnitude Spectrum"),plt.xticks([]),plt.yticks([])
plt.show()

可以看到输出结果的中心部分更白（亮），这说明低频分量更多。
现在可以进行频域变换了，可以在频域对图像进行一些操作了，例如高通滤波和重构图像（DFT的逆变换），比如可以使用一个60x60的矩形窗口对图像进行掩模操作从而去除低频分量。然后再使用函数np.fft.ifftshift()进行逆平移操作，所以现在直流分量有回到左上角了，最后使用函数np.ifft2()进行FFT逆变换。同样又得到一堆复杂的数字，然后取绝对值：

import cv2 
import numpy as np 
from matplotlib import pyplot as plt
img = cv2.imread('image/rose.jpg', 0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
rows, cols = img.shape
crow, ccol = int(rows/2), int(cols/2)
fshift[crow-30:crow+30, ccol-30:ccol+30] = 0
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)
plt.subplot(131),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(132),plt.imshow(img_back, cmap = 'gray')
plt.title('Image after HPF'), plt.xticks([]), plt.yticks([])
plt.subplot(133),plt.imshow(img_back)
plt.title('Result in JET'), plt.xticks([]), plt.yticks([])
plt.show()

结果显示高通滤波其实就是一种边界检测操作。上图的JET会有振铃效应。这是由于使用矩形窗口做掩模造成的。这个掩模被转换成正弦形状就会出现，所有一般不使用矩形窗口滤波，最好的选择是高斯窗口；

2.2 OpenCV中的傅里叶变换
OpenCV中响应的函数是cv2.dft()和cv2.idft()。和前面的输出结果一样。但是爽通道的，第一个结果是实数部分，第二个结果是虚数部分。输入图像首先要转换成np.float32格式。实例如下：

import cv2 
import numpy as np 
from matplotlib import pyplot as plt
img = cv2.imread('image/rose.jpg', 0)
dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)
magnitude_spectrum = 20*np.log(cv2.magnitude(dft_shift[:,:,0], dft_shift[:,:,1]))
plt.subplot(121),plt.imshow(img, cmap='gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap='gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

注意：可以使用函数cv2.cartToPolar()，它会同事返回振幅和相位。
现在我们来做逆DFT。在前面的部分实现了一个HPF(高通滤波)，现在来做LPF(低通滤波)将高频部分去除。其实就是对图像进行模糊操作。首先需要构建一个掩模，与低频区域对应的地方设置为1，与高频区域对应的地方设置为0。

img = cv2.imread('image/rose.jpg', 0)
dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)
rows, cols = img.shape 
crow, ccol = int(rows/2), int(cols/2)
# create a mask first, center square is 1, remaining all zeros
mask = np.zeros((rows, cols, 2), np.uint8)
mask[crow-30:crow+30, ccol-30:ccol+30] = 1
# apply mask and inverse DFT
fshift = dft_shift * mask
f_ishift = np.fft.ifftshift(fshift)
img_back = cv2.idft(f_ishift)
img_back = cv2.magnitude(img_back[:,:,0], img_back[:,:,1])
plt.subplot(121),plt.imshow(img, cmap='gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(img_back, cmap='gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

注意：OpenCV中的函数cv2.dft()和cv2.idft()要比Numpy快；当数组的大小为某些值时DFT的性能会更好。当数组的大小是2的指数时DFT效率最高。当数组的大小是2,3,5的倍数时效率也会很高。所以如果你想要提高代码的运行效率时，可以修改输入图像的大小（补0）。对于OpenCV你必须自己动手补0.但是Numpy，只需要指定FFT运算的大小，它就会自动补0。
那么如何确定最佳大小呢？OpenCV提供了一个函数：cv2.getOptimalDFTSize()。他可以同事被cv2.dft()和np.fft.fft2()使用；

3、 模板匹配
模板匹配是用来在一副大图中搜寻查找模板图像位置的方法。OpenCV为我们提供了函数：cv2.matchTemplate()和2D卷积一样，它也是用图像在输入图像（大图）上滑动，并在每一个位置对模板图像和与其对应的输入图像的子区域进行比较。返回的结果是一个灰度图像，每一个像素值表示了此区域与模板的匹配程度。
如果输入图像的大小是W*H，模板的大小是w*h，输出的结果的大小就是W-w+1*H-h+1；当得到这幅图之后，就可以使用函数cv2.minMaxLoc()来找到其中的最小值和最大值的位置了。第一个值为矩形左上角的点的位置，(w,h)为模板矩形的宽和高，这个矩形就是找到的模板区域了。注意：如果使用的比较方法是cv2.TM_SQDIFF,最小值对应的位置才是匹配的区域。
3.1 OpenCV中的模板匹配

# -*- coding:utf-8 -*- 

import cv2 
import numpy as np 
from matplotlib import pyplot as plt 

img = cv2.imread("./image/nba.jpg", 0)
img2 = img.copy()
template = cv2.imread("./image/weide.jpg", 0)
w, h = template.shape[::-1]

# All the 6 methods for comparison in a list 
methods = ['cv2.TM_CCOEFF', 'cv2.TM_CCOEFF_NORMED', 'cv2.TM_CCORR',
            'cv2.TM_CCORR_NORMED', 'cv2.TM_SQDIFF', 'cv2.TM_SQDIFF_NORMED']

for meth in methods:
    img = img2.copy()
    # exec 语句用来执行存储在字符串或者文件中的Python语句
    # eval 语句用来计算存储在字符串中的有效Python表达式
    method = eval(meth)
    # Apply template Matching 
    res = cv2.matchTemplate(img, template, method)
    min_val, max_val, min_loc, max_loc = cv2.minMaxLoc(res)
    # 使用不同的比较方法，对结果的解释不同
    # if the method is TM_SQDIFF or TM_SQDIFF_NORMED, take minimum 
    if method in [cv2.TM_SQDIFF, cv2.TM_SQDIFF_NORMED]:
        top_left = min_loc
    else:
        top_left = max_loc

    bottom_right = (top_left[0]+w, top_left[1]+h)
    cv2.rectangle(img, top_left, bottom_right, 255 , 2)

    plt.subplot(121),plt.imshow(res, cmap='gray')
    plt.title('Matching Result'),plt.xticks([]),plt.yticks([])
    plt.subplot(122),plt.imshow(img, cmap='gray')
    plt.title('Detected Point'), plt.xticks([]),plt.yticks([])
    plt.suptitle(meth)

    plt.show()

3.2 对象的模板匹配
目标对象在图像中出现很多次怎么办？函数cv2.imMaxLoc()只会给出最大值和最小值。此时就要使用阈值了。

import cv2 
import numpy as np 
from matplotlib import pyplot as plt 

img_rgb = cv2.imread('./image/nba.jpg')
img_gray = cv2.cvtColor(img_rgb, cv2.COLOR_BGR2GRAY)
template = cv2.imread('./image/weide.jpg', 0)
w, h = template.shape[::-1]

res = cv2.matchTemplate(img_gray, template, cv2.TM_CCOEFF_NORMED)
threshold = 0.8 

loc = np.where(res >= threshold)
for pt in zip(*loc[::-1]):
    cv2.rectangle(img_rgb, pt, (pt[0]+w, pt[1]+h), (0,0,255), 2)

cv2.imwrite('res.png', img_rgb)

4. 分水岭算法图像分割
任何一副灰度图像都可以被看成拓扑平面，灰度值高的区域可以被看成是山峰，灰度值低的区域可以被看成是山谷。向每一个山谷中灌不同颜色的水。随着水位的升高不同山谷的水就会相遇汇合，为了防止不同山谷的水汇合，需要在水汇合的地方构建起堤坝，不停的灌水，不停的构建堤坝直到所有的山峰被水淹没。我们构建好的堤坝就是对图像的分割，这就是分水岭算法。http://cmm.ensmp.fr/~beucher/wtshed.html （此链接是动画演示）
但是这种方法通常都会得到过度分割的结果，这是由噪声或者图像中其他不规律的因素造成的。为了减少这种影响，OpenCV采用了基于掩模的分水岭算法，在这中算法中药设置哪些山谷点会汇合，哪些不会。这是一种交互式的图像分割。要做的就是给我们已知的最想打上不同的标签。如果某个区域肯定不是对象而是背景就是用另外一个颜色标签标记。而剩下的不能确定是前景还是背景的区域用0标记，这就是标签。然后实施分水岭算法。每一个灌水，标签就会被更新，当两个不同颜色的标签相遇时就构建堤坝，知道所有的山峰掩模，最后得到的边界对象（堤坝）的值为-1。实例如下：

# -*- coding:utf-8 -*-
import numpy as np 
import cv2 
from matplotlib import pyplot as plt 

img = cv2.imread('./image/bi.jpg')
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
ret, thresh = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY_INV+cv2.THRESH_OTSU)
cv2.imshow('1', thresh)
# noise removal 
kernel = np.ones((3,3), np.uint8)
opening = cv2.morphologyEx(thresh, cv2.MORPH_OPEN, kernel, iterations=2)
# sure background area
sure_bg = cv2.dilate(opening, kernel, iterations=3)
# finding sure foreground area 
# 距离变换的基本含义是计算一个图像中非零像素点到最后的零像素点的距离
# 也就是到零像素点的最短距离
# 最常见的距离变换算法就是通过连续的腐蚀操作来实现，腐蚀操作的停止条件
# 是所有前景像素都被完全腐蚀。这样根据腐蚀的先后顺序，我们就得到各个前景像素
# 点到前景中心固件像素点的距离，根据各个像素点的距离值，设置为不同的灰度值
# 这就完成了二值图像的距离变换
# cv2.distanceTransform(src, distanceType, maskSize)
# 第二个参数0，1，2分别代表CV_DIST_L1,CV_DIST_L2,CV_DIST_C
dist_transform = cv2.distanceTransform(opening,1,5)
ret, sure_fg = cv2.threshold(dist_transform, 0.7*dist_transform.max(), 255,0)
# Finding unknown region
sure_fg = np.uint8(sure_fg)
unknown = cv2.subtract(sure_bg, sure_fg)
# marker labelling 
ret, markers1 = cv2.connectedComponents(sure_fg)
# add one to all labels so that sure background is not 0, but 1
markers = markers1+1
# Now, mark the region of unkown with zero
markers[unknown==255] = 0

markers3 = cv2.watershed(img, markers)
img[markers3 == -1] = [255, 0, 0]

cv2.imshow('zuihou', img)
cv2.waitKey(0)
cv2.destroyAllWindows()

5. 使用GrabCut算法进行交互式前景提取
开始时用户需要用一个矩形将前景区域框住（前景区域应该完全被包括在矩形框内部）。然后算法进行迭代式分割直达到最好的结果。但是有时分割的结果不够理想，比如吧前景当成了背景，或者把背景当成了前景。这种情况下，就需要用户来进行修改了，用户只需要在不理想的部位修改一下就可以了。
过程如下：
1. 用户输入一个矩形，矩形外的所有区域坑定都是背景。矩形内的东西是未知的，同样用户确定前景和背景的任何操作都不会被程序改变。
2. 计算机会对我们的输入图像做一个初始化标记。它会标记前景和背景像素。
3. 使用一个高斯混合模型（GMM）对前景和背景建模。
4. 根据我们的输入，GMM会学习并创建新的像素分布。对那些分类位置的像素（可能是前景也可能是背景），可以根据他们与已知分类（如背景）的像素的关系来进行分类（就像是在做聚类操作）
5. 这样就会根据像素的分布创建一副图。图中的节点就是像素点。除了像素点做节点之外还有两个节点：Source_node和Sink_node。所有的前景像素都和Source_node相连。所有的背景像素和Sink_node相连。
6. 将像素连接到Source_node/end_node的（边）的权重由它们属于同一类（同是前景或同是背景）的概率来决定。两个像素之间的权重由边的信息或者两个像素的相似性来决定。如果两个像素的颜色由很大的不同，那么它们之间的边的权重就会很小。
7. 使用mincut算法对上面得到的图进行分割。它会根据最低成本方程将图分为Source_node和Sink_node。成本方程就是被减掉的所有边的权重之和。在裁剪之后，所有连接到Source_node的像素被认为是前景，所有连接到Sink_node的像素被认为是背景。
8. 继续这个过程知道分类收敛。

OpenCV提供了函数：cv2.grabCut()；先看看参数：
1. img-输入图像
2.mask-掩模图像，用来确定哪些区域是背景，前景，可能是前景/背景等。可以设置为cv2.GC_BGD,cv2.GC_FGD,cv2.GC_PR_BGD,cv2.GC_PR_FGD,或者直接输入0,1，2,3也行。
3. rect-包含前景的矩形。格式为（x,y,w,h）
4. bdgModel,fgdModel-算法内部使用的数组，你只需要创建两个大小为(1,65)，数据类型为np.float64的数组。
5. iterCount-算法的迭代次数
6. mode可以设置为cv2.GC_INIT_WITH_RECT或cv2.GC_INIT_WITH_MASK，也可以联合使用。这是用来确定我们进行修改的方式，矩形模式或者掩模模式
import numpy as np
import cv2
from matplotlib import pyplot as plt

img = cv2.imread('./image/messi5.jpg')
mask = np.zeros(img.shape[:2], np.uint8)
bgdModel = np.zeros((1,65), np.float64)
fgdModel = np.zeros((1,65), np.float64)

rect = (50, 50, 450, 290)
# 函数的返回值是更新的mask,bgdModel,fgdModel
cv2.grabCut(img, mask, rect, bgdModel, fgdModel, 5, cv2.GC_INIT_WITH_RECT)

newmask = cv2.imread('./image/newmessi5.jpg', 0)
mask[newmask == 0] = 0
mask[newmask == 255] = 1
mask, bgdModel, fgdModel = cv2.grabCut(img,mask,None,bgdModel,fgdModel,5,cv2.GC_INIT_WITH_MASK)

mask = np.where((mask==2)|(mask==0), 0, 1).astype('uint8')
img = img*mask[:,:,np.newaxis]
b, g, r = cv2.split(img)
img = cv2.merge([r,g,b])
plt.imshow(img)
plt.colorbar()
plt.show()

结语：本章节有点多，但都是很必要的知识点，下一章节会讲图像的直方图，变换，模板匹配以及Hough变换和图像分割。谢谢