数字图像处理

主要参考自：OpenCV中文官方文档。

一、OpenCV相关

1. Opencv中的坐标系与Numpy矩阵坐标系

在OpenCV中坐标系横向为x轴，纵向为y轴；而在Numpy矩阵坐标系中，纵向为x轴，横向为y轴。因此，在OpenCV中设置处理的相关坐标时，例如绘制图形、Resize时，应该是采用(w,h)进行设置，但是在访问image中的数值时或者是获取shape等Numpy操作时，此时image为Numpy矩阵，应该使用(h,w)进行访问。
更详细的分析参见：OpenCV 中图像坐标系统与Python中NumPy Arrays之间的关系

2. 阈值化

（1）简单阈值

（2）自适应阈值

在上一节中，我们使用一个全局值作为阈值。但这可能并非在所有情况下都很好，例如，如果图像在不同区域具有不同的光照条件。在这种情况下，自适应阈值阈值化可以提供帮助。在此，算法基于像素周围的小区域确定像素的阈值。因此，对于同一图像的不同区域，我们获得了不同的阈值，这为光照度变化的图像提供了更好的结果。

除上述参数外，方法cv.adaptiveThreshold还包含三个输入参数：
该adaptiveMethod决定阈值是如何计算的：
cv.ADAPTIVE_THRESH_MEAN_C::阈值是邻近区域的平均值减去常数C。 cv.ADAPTIVE_THRESH_GAUSSIAN_C:阈值是邻域值的高斯加权总和减去常数C。
该BLOCKSIZE确定附近区域的大小，C是从邻域像素的平均或加权总和中减去的一个常数。

下面的代码比较了光照变化的图像的全局阈值和自适应阈值：

（3）Otsu的二值化

在全局阈值化中，我们使用任意选择的值作为阈值。相反，Otsu的方法避免了必须选择一个值并自动确定它的情况。

考虑仅具有两个不同图像值的图像（双峰图像），其中直方图将仅包含两个峰。一个好的阈值应该在这两个值的中间。类似地，Otsu的方法从图像直方图中确定最佳全局阈值。

为此，使用了cv.threshold作为附加标志传递。阈值可以任意选择。然后，算法找到最佳阈值，该阈值作为第一输出返回。

查看以下示例。输入图像为噪点图像。在第一种情况下，采用值为127的全局阈值。在第二种情况下，直接采用Otsu阈值法。在第三种情况下，首先使用5x5高斯核对图像进行滤波以去除噪声，然后应用Otsu阈值处理。了解噪声滤波如何改善结果。

import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt
img = cv.imread('noisy2.png',0)
# 全局阈值
ret1,th1 = cv.threshold(img,127,255,cv.THRESH_BINARY)
# Otsu阈值
ret2,th2 = cv.threshold(img,0,255,cv.THRESH_BINARY+cv.THRESH_OTSU)
# 高斯滤波后再采用Otsu阈值
blur = cv.GaussianBlur(img,(5,5),0)
ret3,th3 = cv.threshold(blur,0,255,cv.THRESH_BINARY+cv.THRESH_OTSU)
# 绘制所有图像及其直方图
images = [img, 0, th1,
          img, 0, th2,
          blur, 0, th3]
titles = ['Original Noisy Image','Histogram','Global Thresholding (v=127)',
          'Original Noisy Image','Histogram',"Otsu's Thresholding",
          'Gaussian filtered Image','Histogram',"Otsu's Thresholding"]
for i in xrange(3):
    plt.subplot(3,3,i*3+1),plt.imshow(images[i*3],'gray')
    plt.title(titles[i*3]), plt.xticks([]), plt.yticks([])
    plt.subplot(3,3,i*3+2),plt.hist(images[i*3].ravel(),256)
    plt.title(titles[i*3+1]), plt.xticks([]), plt.yticks([])
    plt.subplot(3,3,i*3+3),plt.imshow(images[i*3+2],'gray')
    plt.title(titles[i*3+2]), plt.xticks([]), plt.yticks([])
plt.show()

Otsu的阈值的查找

二、数字图像处理基础

1. 仿射变换

仿射变换矩阵的形式：

常见变换的矩阵：

以上图片来自于：仿射变换及其变换矩阵的理解

2. 滤波

低通滤波器（LPF），高通滤波器（HPF）等对图像进行滤波。LPF有助于消除噪声，使图像模糊等。HPF滤波器有助于在图像中找到边缘。

（1）均值滤波

通过将图像与归一化框滤镜进行卷积来完成的。它仅获取内核区域下所有像素的平均值，并替换中心元素。
3x3归一化框式过滤器如下所示：
$$K=\frac{1}{9}\left[\begin{array}{lll}1& 1& 1\\ 1& 1& 1\\ 1& 1& 1\end{array}\right]$$

（2）高斯模糊

采用像素周围的邻域并找到其高斯加权平均值作为中心元素的值。高斯模糊对于从图像中去除高斯噪声非常有效。

（3）中值模糊

提取内核区域下所有像素的中值，并将中心元素替换为该中值。这对于消除图像中的椒盐噪声非常有效。在中值模糊中，中心元素总是被图像中的某些像素值代替。有效降低噪音。其内核大小应为正奇数整数。

（4）双边滤波

在去除噪声的同时保持边缘清晰锐利非常有效。高斯滤波器采用像素周围的邻域并找到其高斯加权平均值。高斯滤波器仅是空间的函数，也就是说，滤波时会考虑附近的像素。它不考虑像素是否具有几乎相同的强度。它不考虑像素是否是边缘像素。因此它也模糊了边缘，这是我们不想做的。

双边滤波器在空间中也采用高斯滤波器，但是又有一个高斯滤波器，它是像素差的函数。空间的高斯函数确保仅考虑附近像素的模糊，而强度差的高斯函数确保仅考虑强度与中心像素相似的那些像素的模糊。由于边缘的像素强度变化较大，因此可以保留边缘。

以上截图来自于：双边滤波（Bilateral Filter）详解。

3. 形态学操作

（1）腐蚀操作

它侵蚀前景物体的边界(尽量使前景保持白色)。原始图像中的一个像素(无论是1还是0)只有当内核下的所有像素都是1时才被认为是1，否则它就会被侵蚀(变成0)。

结果是，根据内核的大小，边界附近的所有像素都会被丢弃。因此，前景物体的厚度或大小减小，或只是图像中的白色区域减小。它有助于去除小的白色噪声，分离两个连接的对象等。

（2）膨胀操作

它与腐蚀操作正好相反。如果内核下的至少一个像素为“ 1”，则像素元素为“ 1”。因此，它会增加图像中的白色区域或增加前景对象的大小。通常，在消除噪音的情况下，腐蚀后会膨胀。因为腐蚀会消除白噪声，但也会缩小物体。因此，我们对其进行了扩展。由于噪音消失了，它们不会回来，但是我们的目标区域增加了。在连接对象的损坏部分时也很有用。

（3）开运算

开运算只是腐蚀然后膨胀的另一个名称。如上文所述，它对于消除噪音很有用。

（4）闭运算

闭运算与开运算相反，先膨胀然后再腐蚀。在关闭前景对象内部的小孔或对象上的小黑点时很有用。

（5） 形态学梯度

图像膨胀和腐蚀之间的区别。结果将看起来像对象的轮廓。

（6）白帽

它是输入图像和图像开运算之差。

（7）黑帽

它是输入图像和图像闭运算之差。

4. 图像梯度

常见的梯度滤波器或高通滤波器，有Sobel，Scharr和Laplacian。

（1）Sobel、Scharr算子

Sobel算子垂直方向和水平方向的模板如上图所示，前者可以检测出图像中的水平方向的边缘，后者则可以检测图像中垂直方向的边缘。实际应用中，每个像素点取两个模板卷积的最大值作为该像素点的输出值，运算结果是一副边缘图像。

当kernel_size=3时，使用3x3 Scharr滤波器，比3x3 Sobel滤波器具有更好的结果。

以上参考自索贝尔算子、sobel算子原理与实现、边缘检测 - Scharr 滤波器。

（2）Laplacian 算子

它计算了由关系$\Delta src=\frac{{\partial }^{2}src}{\partial x^2}+\frac{{\partial }^{2}src}{\partial y^2}$给出的图像的拉普拉斯图,它是每一阶导数都是通过Sobel算子计算。如果ksize = 1,然后使用以下内核用于过滤:
$$\text{ kernel }=\left[\begin{array}{ccc}0& 1& 0\\ 1& -4& 1\\ 0& 1& 0\end{array}\right]$$

参考自：为什么 空间二阶导（拉普拉斯算子）这么重要？ - ChainingBlocks的回答 - 知乎。

5. Canny边缘检测

该部分参考自：Canny边缘检测算法、Python - Opencv 之 Canny 边缘检测。

Canny边缘主要包含四个步骤：

（1）高斯滤波

采用高斯滤波器去除部分噪声。

（2）梯度计算

基于Sobel算子计算x以及y方向的梯度，并且计算出各个位置的梯度大小$G$和方向$\theta$：
$$G=\sqrt{G_x^2+G_y^2}$$$$\theta ={\tan }^{-1}\left(\frac{G_y}{G_x}\right)$$

（3）非极大值抑制（Non-Maxim Suppression，NMS）

计算得到梯度值和梯度方向后，对图片进行全面的扫描，以去除不构成边缘的无关像素点。
对于每个像素，检查其是否是在梯度方向中其临近像素点中的局部最大值。如果一个像素点属于边缘，那么这个像素点在梯度方向上的梯度值是最大的。否则不是边缘，将灰度值设为0。NMS会使边缘尽可能薄。

（4）阀值过滤

需要设定两个阈值，minVal 和 maxVal。
任何边缘的强度梯度大于 maxVal 的确定为边缘；而小于 minVal 的确定为非边缘，并丢弃。
位于 maxVal 和 minVal 阈值间的边缘为待分类边缘，或非边缘，此时再基于连续性进行判断。如果边缘像素连接着 “确定边缘” 像素，则认为该边缘属于真正边缘的一部分；否则，丢弃该边缘。

6. 图像金字塔

参考自：图像金字塔（高斯金字塔、拉普拉斯金字塔）。

一幅图像的金字塔是一系列以金字塔形状排列的分辨率逐步降低，且来源于同一张原始图的图像集合。层级越高，则图像越小，分辨率越低。
主要有两种类型的图像金字塔：

• 高斯金字塔：用于向下采样。
• 拉普拉斯金字塔：用于从金字塔高层（分辨率低）重建下层未采样的图像，一般配合高斯金字塔使用。

（1）高斯金字塔的构建

1. 首先使用高斯核对原图进行滤波操作。
2. 然后将所有的偶数行和偶数列去除，则图像的尺寸就会下降为原图的1/4。

反复迭代这个过程就能形成不同尺度的图像，从而形成高斯金字塔。

（2）拉普拉斯金字塔的构建

1. 首先在图像的各个方向上进行扩充，新增的行和列以0进行填充。
2. 然后用构建高斯金字塔同样的内核（乘以四）与放大后的图像进行卷积，从而为新增元素添加近似值，得到近似的放大后的图像。
3. 最后用于放大后的图像一个尺度的高斯金字塔图像减去近似的放大后的图像，即可得到拉普拉斯金字塔当前尺度的图像。

反复迭代这个过程就能形成不同尺度的图像，从而形成拉普拉斯金字塔。
整个拉普拉斯金字塔的运算过程如下图所示：

金字塔的一种应用是图像融合。具体可参见：使用金字塔进行图像融合。

7. 直方图

主要参考自：直方图均衡化、直方图-2：直方图均衡。

（1）基础知识

直方图的横坐标表示一系列的范围，不一定总是0-255的范围，纵坐标是在这些范围内像素点的数目。
这只是理解图像的另一种方式。通过查看图像的直方图，可以直观地了解该图像的对比度，亮度，强度分布等。

• BINS
  对于直方图来说，其最终要的一个参数就是BINS，即横坐标的个数，实际上就是指定范围所需要的划分的范围的个数。

  上面的直方图显示每个像素值的像素数，即从0到255。即，横坐标需要256个值来显示上面的直方图。如果不需要分别找到所有像素值的像素数，而是找到像素值间隔中的像素数怎么办？ 例如，需要找到介于0到15之间的像素数，然后找到16到31之间，…，240到255之间的像素数。只需要16个值即可表示直方图。

  因此，我们要做的就是将整个直方图分成16个子部分，每个子部分的值就是其中所有像素数的总和。 每个子部分都称为“ BIN”。在第一种情况下，bin的数量为256个（每个像素一个），而在第二种情况下，bin的数量仅为16个。BINS由OpenCV文档中的histSize术语表示。

• DIMS
  需要计算直方图的某个维度。

• RANAGE
  需要划分的指定的范围的大小。

（2）计算方法

a. OpenCV中的直方图计算

cv.calcHist（images，channels，mask，histSize，ranges [，hist [，accumulate]]）

• images：它是uint8或float32类型的源图像。它应该放在方括号中，即“ [img]”。
• channels：也以方括号给出。它是我们计算直方图的通道的索引。例如，如果输入为灰度图像，则其值为[0]。对于彩色图像，您可以传递[0]，[1]或[2]分别计算蓝色，绿色或红色通道的直方图。
• mask：图像掩码。为了找到完整图像的直方图，将其指定为“无”。但是，如果要查找图像特定区域的直方图，则必须为此创建一个掩码图像并将其作为掩码。
• histSize：这表示我们的BIN计数。需要放在方括号中。对于全尺寸，我们通过[256]。
• ranges：这是我们的RANGE。通常为[0,256]。

img = cv.imread('home.jpg',0)
hist = cv.calcHist([img],[0],None,[256],[0,256])

hist是256x1的数组，每个值对应于该图像中具有相应像素值的像素数。

b. numpy的直方图计算

Numpy还提供了一个函数np.histogram()。

hist,bins = np.histogram(img.ravel(),256,[0,256])

hist与我们之前计算的相同。但是bin将具有257个元素，因为Numpy计算出bin的范围为0-0.99、1-1.99、2-2.99等。因此最终范围为255-255.99。为了表示这一点，他们还在最后添加了256。但是我们不需要256。最多255就足够了。

另外 Numpy还有另一个函数np.bincount()，它比np.histogram()快10倍左右。因此，对于一维直方图，您可以更好地尝试一下。不要忘记在np.bincount中设置minlength = 256。例如，hist = np.bincount(img.ravel()，minlength = 256)
注意 OpenCV函数比np.histogram()快大约40倍。因此，尽可能使用OpenCV函数。

（3）直方图均衡化

考虑这样一个图像，它的像素值仅局限于某个特定的值范围。例如，较亮的图像将把所有像素限制在高值上。但是一幅好的图像会有来自图像所有区域的像素。因此，您需要将这个直方图拉伸到两端(如下图所示，来自wikipedia)，这就是直方图均衡化的作用(简单来说)。这通常会提高图像的对比度。

a. 简单的直方图均衡化

简单总结一下直方图均衡化的步骤：

• 计算直方图的累积分布
• 根据累积分布的最小以及最大值，对其进行归一化，并映射至[0, 255]的范围，该映射即为均衡化的映射函数
• 将上述映射函数应用于原图

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

img = cv.imread('wiki.jpg',0)
hist,bins = np.histogram(img.flatten(),256,[0,256])  # 计算直方图
cdf = hist.cumsum()  # 计算直方图的累积分布

# 去除掉出现0次的分布，然后进行最大最小值归一化，建立映射
cdf_m = np.ma.masked_equal(cdf,0)
cdf_m = (cdf_m - cdf_m.min())*255/(cdf_m.max()-cdf_m.min())
cdf = np.ma.filled(cdf_m,0).astype('uint8')

# 将映射应用于原图
# 以下代码的意义是：将img的某个位置的对应的值作为索引，去查询cdf对应的值对应的映射值，将该值作为该位置的新值，从而形成img2
img2 = cdf[img] 

b. CLAHE（对比度受限的自适应直方图均衡）

简单的直方图均衡化考虑了图像的整体对比度，但这样的一个坏处在于可能会降低重要区域的对比度。因此，提出了自适应直方图均衡化，即在一个个的小块中进行直方图均衡化。
在较小的区域中，直方图将限制在一个较小的区域中（除非存在噪声）。如果有噪音，它将被放大。为了避免这种情况，应用了对比度限制。如果任何直方图bin超出指定的对比度限制（在OpenCV中默认为40），则在应用直方图均衡之前，将这些像素裁剪并均匀地分布到其他bin。

8. 霍夫变换

（1）霍夫线变换

一条线可以表示为$y=mx+c$或以参数形式表示为$\rho =x\cos \theta +y\sin \theta$，其中$\rho$是从原点到该线的垂直距离，而$\theta$是由该垂直线和水平轴形成的角度以逆时针方向测量。

如果线在原点下方通过，则它将具有正的$\rho$且角度小于180。如果线在原点上方，则将角度取为小于180，而不是大于180的角度。$\rho$取负值。任何垂直线将具有0度，水平线将具有90度。

现在，让我们看一下霍夫变换如何处理线条。任何一条线都可以用$(\rho ,\theta )$这两个参数表示。
因此，首先创建2D数组或累加器（以保存两个参数的值），并将其初始设置为0。行表示$\rho$，列表示$\theta$。
阵列的大小取决于所需的精度。假设您希望角度的精度为1度，则需要180列。对于$\rho$，最大距离可能是图像的对角线长度。因此，以一个像素精度为准，行数可以是图像的对角线长度。

• 对于待检测的图像，其首先一般需要经过Canny滤波等方式提取边界，降低噪声，得到一个二值化的图像。
• 对于该二值化图像上的每个不为0的像素点，即所有前景像素，进行遍历，对于每个前景像素来说，经过该位置$(x_0,y_0)$的所有的直线的方程可以表示为一簇函数：
  $$\rho =x_0\cdot \cos (\theta )+y_0\cdot \sin (\theta )$$
• 该簇函数在$(\rho ,\theta )$阵列上所有经过的位置对应的累加器的数值都加1。
• 对于所有前景像素点完成上述遍历，并更新累加器的值，该累加器的值实际上意味着该累加器对应的直线在待检测的图像上的像素数目。
• 最后通过一个阈值，输出累加器数目大于一定阈值的直线，即为最终检测到的直线。

OpenCV相关代码为：

import cv2 as cv
import numpy as np
img = cv.imread(cv.samples.findFile('sudoku.png'))
gray = cv.cvtColor(img,cv.COLOR_BGR2GRAY)
edges = cv.Canny(gray,50,150,apertureSize = 3)
lines = cv.HoughLines(edges,1,np.pi/180,200)
for line in lines:
    rho,theta = line[0]
    a = np.cos(theta)
    b = np.sin(theta)
    x0 = a*rho
    y0 = b*rho
    x1 = int(x0 + 1000*(-b))
    y1 = int(y0 + 1000*(a))
    x2 = int(x0 - 1000*(-b))
    y2 = int(y0 - 1000*(a))
    cv.line(img,(x1,y1),(x2,y2),(0,0,255),2)
cv.imwrite('houghlines3.jpg',img)

HoughLines参数说明：

1. 第一个参数，输入图像应该是二进制图像，因此在应用霍夫变换之前，请应用阈值或使用Canny边缘检测。
2. 第二和第三参数分别是ρ和θ精度。
3. 第四个参数是阈值，这意味着应该将其视为直线的最低投票。它表示应检测到的最小线长。

参考资料：OpenCV图像处理|1.20 霍夫线变换、霍夫线变换。

（2）霍夫圆变换

霍夫圆变换与霍夫直线变换原理类似，每个圆对应的$(x-a{)}^2+(y-b{)}^2=r^2$方程可以转化为
$$\begin{array}{l}a=x-r\cdot \cos \theta \\ b=y-r\cdot \sin \theta \end{array}$$
也就是说，此时，我们需要$(x,y,r)$三个参数去确定一个圆。

OpenCV中采用的“霍夫梯度法”以降低复杂度，具体参见：霍夫圆变换原理