数字图像处理笔记

图像的存储

图像有单通道和多通道之分，图像按照多维数组方式存储，如这是一种图像存储方式，Width为单通道的宽，height为单通道的高，channel为通道数

image[nWidth][nHeight][nChannel];

opencv 图像存储格式

opencv中，存储图像使用的是Mat结构（也是一个矩阵）。 初始化可以用：

cv::Mat mat;
cv::Mat mat(rows, cols, type);
cv::Mat mat = cv::Mat::zeros(rows, cols, type);

type: 数据类型和通道数，常用类型如下：

• CV_8UC1：单通道 8 位无符号整数（灰度图）。
• CV_8UC3：三通道 8 位无符号整数（彩色图像，如 BGR）。
• CV_32FC1：单通道 32 位浮点数（常用于计算机视觉中的图像处理）。

Mat的主要成员如下：

• data:指向矩阵数据的指针。
• rows/cols/channels(): 行数、列数、通道数
• type():返回矩阵的数据类型和通道数
• clone()、copyTo(newmat):数据全复制

仿射变换

仿射变换 = 线性变换（旋转、切边、缩放） + 平移

• 缩放 \(\left( \begin{matrix} x'\\ y' \end{matrix} \right) = \left( \begin{matrix} s_x & 0\\ 0 & s_y \end{matrix} \right) \left( \begin{matrix} x\\ y \end{matrix} \right)\)

• 切变 \(\left( \begin{matrix} x'\\ y' \end{matrix} \right) = \left( \begin{matrix} 1 & a\\ 0 & 1 \end{matrix} \right) \left( \begin{matrix} x\\ y \end{matrix} \right)\)

• 旋转 \(\left( \begin{matrix} x'\\ y' \end{matrix} \right) = \left( \begin{matrix} cos\theta & -sin\theta\\ sin\theta & cos\theta \end{matrix} \right) \left( \begin{matrix} x\\ y \end{matrix} \right)\)

空间域图像处理

空间域：图像平面本身，图像的每个像素单元

空间域图像处理包括：灰度变换、空间域滤波

灰度变换

常见的灰度变换：

• 图像灰度反转 \(s = L - 1 - r\)
• 对数变换 \(s = c\cdot log_{v+1}(1+v\cdot r)\)
• Gamma变换 \(s = cr^\gamma\)
• 直方图均衡化

空间域滤波

操作原理：给定待处理图像、核（kernel）。将核与待处理图像对应位置相乘再求平均得到的数值即是得到的结果（卷积）

平滑

平滑滤波器包括：

• 平滑线性滤波器（盒子滤波器）

所有系数都相等。核矩阵的尺寸越大，滤波结果越模糊

opencv库函数为:blur,使用方法：

for i in range(1, MAX_KERNEL_LENGTH, 2):
    dst_img = cv2.blur(img, (i, i))

• 中值滤波器

像素邻域内灰度的中值替换目标像素，用于去除椒盐噪声

中值滤波器没有固定的核矩阵结构，是将邻域中的所有数取中位数得到结果

opencv库函数为:medianBlur,使用方法：

for i in range(1,MAX_KERNEL_LENGTH , 2):
    dst_img = cv2.medianBlur(img, i)

• 高斯滤波器

是一种线性滤波器，有效抑制噪声，平滑图像。计算方式与平滑线性滤波器相同，区别为核结构不同。

对于大小为(2k+1) x (2k+1)的高斯滤波器，核矩阵的各个元素按如下公式给出：

标准差$\sigma$越大，平滑效果越明显，越小越比明显，取自身的权重就越大 opencv库函数为:GaussianBlur,使用方法：

for i in range(1, MAX_KERNEL_LENGTH, 2):
    dst_img = cv2.GaussianBlur(img, (i, i), 0)

• 双边滤波器

是一种非线性滤波器，基本思想是同时考虑空间邻域和 像素值差异两个因素来平滑图像。能够在去噪的同时保留图像的边缘信息。

opencv库函数为:bilateralFilter,使用方法：

for i in range(1, MAX_KERNEL_LENGTH, 2):
    dst_img = cv2.GaussianBlur(img, i, i * 2, i / 2)

空间域平滑滤波作用：

1. 模糊处理：大目标提取之前先去掉图像中的琐碎细节
2. 降低噪声：去掉灰度的急剧变化

锐化

作用：突出灰度的过度部分（寻找边界、增强图像）

原理：对图像求一阶导数和二阶导数

锐化滤波器包括：

• Prewitt算子：可以准确定位到不同类型的边缘线

• Sobel算子：比Prewitt算子更突出中心点达到平滑目的

• Laplace算子：可以检测到灰度突变的点

将原图像加上或减去锐化图像结果可以得到增强后的图像

opencv中图像锐化写法：

img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
kernel = np.array([[0,-1,0],[-1,4,-1],[0,-1,0]], dtype=int)
dst_img = cv2.filter2D(img_gray, -1, kernel)

频率域滤波

步骤：

• 给定一幅大小为$M\times N$的输入图像，从$P>=2M-1$和$Q>=2N-1$得到填充参数P和Q，典型地，选择P=2M和Q=2N。
• 对$f(x,y)$添加必要数量的0，形成大小为$P\times Q$的
• 用$(-1)^{x+y}$乘以$f_p(x,y)$移到其变换的中心
• 计算来自步骤3的图像的DFT，得到$F(u,v)$
• 生成一个真实的、对称的滤波函数$H(u,v)$,其大小为$P \times Q$，中心在$(p/2,Q/2)$处，用阵列相乘形成乘积$G(u,v)=H(u,v)F(u,v)$；即$G(i,k)=H(i,k)F(i,k)$
• 得到处理后的图像：$g_p(x,y)={real[F’^{-1}[G(u,v)]]}(-1)^{x+y}$，其中，为忽略由于计算不准确导致的寄生复分量，选择了实部，下标p指出我们处理的式填充后的序列
• 通过$g_p(x,y)$的左上象限提取$M\times N$区域，得到最终处理结果g(x,y)

opencv中示例代码为：

# 读取图像
img_o = cv2.imread('cat.png')
img = cv2.cvtColor(img_o, cv2.COLOR_BGR2GRAY)

# 获取图像的尺寸
M, N = img.shape

# 设置 P 和 Q 为图像尺寸的两倍
P = 2 * M
Q = 2 * N

# 步骤 1: 填充图像到 P x Q
# 在右边和下边填充零，使图像大小变为 P x Q
f_p = np.zeros((P, Q), dtype=np.float32)
f_p[:M, :N] = img

# 步骤 2: 创建 (-1)^(x + y) 的矩阵
f_p_shifted = np.copy(f_p)  # 创建一个与 f_p 相同的矩阵

for i in range(P):  # 遍历矩阵的每个元素，进行相乘
    for j in range(Q):
        f_p_shifted[i, j] *= (-1) ** (i + j)

# 步骤 3: 对图像进行傅里叶变换
F_uv = np.fft.fft2(f_p_shifted)

# 步骤 4: 生成一个对称的滤波函数 H(u, v),这里使用简单的理想低通滤波器为例
u, v = np.meshgrid(np.fft.fftfreq(P), np.fft.fftfreq(Q))
distance = np.sqrt(u ** 2 + v ** 2)
H_uv = np.where(distance < 0.7, 1, 0)  # 设置滤波器大小，低通部分为1，高通部分为0

# 步骤 5: 滤波操作
G_uv = H_uv * F_uv

# 步骤 6: 逆傅里叶变换
g_p_shifted = np.fft.ifft2(G_uv)

# 反向移动频谱中心
g_p = np.real(g_p_shifted) * (-1) ** (np.indices(g_p_shifted.shape).sum(axis=0))

# 步骤 7: 提取左上象限，恢复原图大小
g = g_p[:M, :N]

平滑

理想低通滤波器

会产生振铃效应，一般不用

u, v = np.meshgrid(np.fft.fftfreq(P), np.fft.fftfreq(Q))
distance = np.sqrt(u ** 2 + v ** 2)
H_uv = np.where(distance < D_0, 1, 0)  # 设置滤波器大小，低通部分为1，高通部分为0

巴特沃斯Butterworth低通滤波器

此滤波器可以严格控制截止频率在低频和高频之间的过渡，产生轻微振铃效应

n为滤波器阶数，阶数越高振铃效应越明显。D_n为截止频率。

u, v = np.meshgrid(np.fft.fftfreq(P), np.fft.fftfreq(Q))
D_uv = np.sqrt((u - P/2) ** 2 + (v - Q/2) ** 2)
H_uv = 1 / (1 + (D_uv / D_n) ** (2 * n))  

高斯Gaussion低通滤波器

没有振铃效应，平滑效果比二阶BLPF差

u, v = np.meshgrid(np.fft.fftfreq(P), np.fft.fftfreq(Q)) 
D_uv = np.sqrt((u - P/2) ** 2 + (v - Q/2) ** 2) 
H_uv = np.exp(-(D_uv ** 2) / (2 * (D_n ** 2)))

锐化

理想高通滤波器

写法上与理想低通滤波器相反，但也有振铃效果

u, v = np.meshgrid(np.fft.fftfreq(P), np.fft.fftfreq(Q))
distance = np.sqrt(u ** 2 + v ** 2)
H_uv = np.where(distance < D_0, 0, 1)  # 设置滤波器大小，低通部分为1，高通部分为0

巴特沃斯Butterworth高通滤波器

此滤波器可以严格控制截止频率在低频和高频之间的过渡，产生轻微振铃效应

n为滤波器阶数，阶数越高振铃效应越明显。D_n为截止频率。

u, v = np.meshgrid(np.fft.fftfreq(P), np.fft.fftfreq(Q))
D_uv = np.sqrt((u - P/2) ** 2 + (v - Q/2) ** 2)
H_uv = 1 - 1 / (1 + (D_0 / D_uv) ** (2 * n))  

高斯Gaussion低通滤波器

没有振铃效应，锐化效果比二阶BLPF差

u, v = np.meshgrid(np.fft.fftfreq(P), np.fft.fftfreq(Q)) 
D_uv = np.sqrt(u ** 2 + v ** 2) 
H_uv = 1 - np.exp(-(D_uv ** 2) / (2 * (D_0 ** 2)))

形态学处理

处理原理： 首先对原图边界进行扩充，使结构元B完全包括在A中（即当结构元原点与原图的每个像素重合时都在新图像中）。

通过结构元对原图像进行集合运算

结构元：研究一幅图像中感兴趣特征所用的小集合或子图像。结构元存在一个原点。

存在不同灰度级的结构元称为不平坦结构元。灰度一致为平坦结构元

腐蚀

对于平坦结构元：将结构元的原点放在图像的每一个像素的位置，对于结构元与图像重合的部分，当结构元所在位置的图像全为1时（子集），该位置标记为1，否则全为0。

公式表示：$A\ominus B = {z|(B)_z\subseteq A}$

opencv函数使用：

erode_img = cv2.erode(img, kernel, iteration = 1) # iteration为腐蚀次数

对于非平坦结构元：对于每一个像素点求结构元与图像像素重合中的每一个像素的差值（图像减结构元）的最小值为该像素的结果。

代码实现：

for i in range(image.shape[0]):
    for j in range(image.shape[1]):
        region = padded_img[i:i + k_h, j:j + k_w]  # 取(扩充过的)图像局部区域
        diff = region - kernel  # 计算像素差值
        erode_img[i, j] = np.min(diff)  # 取最小值作为腐蚀值

膨胀

如果说腐蚀是取“且”，那膨胀是取“或”，但是不同的是结构元需要翻转（按原点对称）

对于平坦结构元：把结构元原点放在图像每个像素的位置，在任何位置的膨胀为包含与结构元翻转（中心对称）重合区域的图像的所有值的最大值

公式表示： $f\oplus b=\mathop{max}\limits_{(s,t)\in \hat b_N} {f(x-s,y-t)}$

opencv函数使用：

dilate_img = cv2.dilate(img, kernel, iteration = 1) # iteration为膨胀次数

作用：填补目标区域中某些空洞以及消除包含在目标区域中的小颗粒噪声

对于非平坦结构元：把结构元原点放在图像每个像素的位置，在任何位置的膨胀为包含与结构元翻转（中心对称）重合区域的图像加上结构元翻转的值的最大值

for i in range(image.shape[0]):
    for j in range(image.shape[1]):
        region = padded_img[i:i + k_h, j:j + k_w]  # 取(扩充过的)图像局部区域
        diff = region + kernel  # 计算像素和值
        dilate_img[i, j] = np.max(diff) # 取最大值作为膨胀值

开运算

先腐蚀后膨胀，腐蚀可以断开微弱连接，然后再膨胀平滑边

公式表示：$f\circ b = (f\ominus b) \oplus b$

opening = cv2.morphologyEx(img, cv2.MORPH_OPEN, kernel)

闭运算

先膨胀，膨胀可以将距离比较近的区域合并起来，再腐蚀掉比较小的区域

公式表示：$f\bullet b = (f\oplus b) \ominus b$

closing = cv2.morphologyEx(img, cv2.MORPH_CLOSE, kernel)

形态学平滑

先进行开运算，再执行闭运算

opening = cv2.morphologyEx(img, cv2.MORPH_OPEN, kernel)
smoothed = cv2.morphologyEx(opening, cv2.MORPH_CLOSE, kernel)

顶帽运算

去掉开操作的结果，即用原图像减去开运算的部分

top_hat = cv2.morphologyEx(img, cv2.MORPH_TOPHAT, kernel)

黑帽运算

闭操作减去原图像

black_hat = cv2.morphologyEx(img, cv2.MORPH_BLACKHAT, kernel)

击中击不中

用一个小的结构元去射击原始图像，击中的元素保留，再用一个很大的结构元去射击图像，击不中的保留，满足击中元素能击中且击不中的元素不能击中的位置就是最终图像

形态学梯度

膨胀结果减去腐蚀结果。强调了区间的边界，同质区域基本不受影响

gradient = cv2.morphologyEx(img, cv2.MORPH_GRADIENT, kernel)

检测与连接

以灰度局部剧烈变化为基础，检测图像中的点、线和边缘，并将检测到的断裂边缘连接起来

检测方法：一阶导数和二阶导数

• 一阶导数产生较粗的边缘
• 二阶导数对于精细细节（如细线、孤立点、噪声）有较强响应
• 二阶导数在灰度斜坡和灰度台阶处会产生双边缘效应
• 二阶导数的符号可以确定是从亮到暗还是从暗到亮

检测孤立点

用二阶导数检测孤立点（拉普拉斯算子）

• 对输入图像，用如下拉普拉斯核进行卷积。

• 对卷积结果进行过滤（阈值处理）：$g(x,y) = \begin{cases} 1\quad if|Z(x,y)|>T
  0\quad otherwise \end{cases}$

检测线

可以用拉普拉斯算子检测，但要根据待检测的线宽设计拉普拉斯核尺寸。对细线和孤立点检测效果较好。但因拉普拉斯核是各向同性的，因此没法检测特定方向线。且对噪声敏感，有双倍加强作用。常产生双像素的边缘

对于不同方向的线，可以用如下的二阶偏导核：

• Horizontal \(\left( \begin{matrix} -1 & -1 & -1\\ 2 & 2 & 2\\ -1 & -1 & -1\\ \end{matrix} \right)\)
• +45°: \(\left( \begin{matrix} 2 & -1 & -1\\ -1 & 2 & -1\\ -1 & -1 & 2\\ \end{matrix} \right)\)
• Vertical: \(\left( \begin{matrix} -1 & 2 & -1\\ -1 & 2 & -1\\ -1 & 2 & -1\\ \end{matrix} \right)\)
• -45°: \(\left( \begin{matrix} -1 & -1 & 2\\ -1 & 2 & -1\\ 2 & -1 & -1\\ \end{matrix} \right)\)

边缘检测

边缘检测基于灰度突变来定位边缘位置。步骤如下：

• 降噪，通过通过高斯滤波器、均值滤波器、中值滤波器等算子对图像进行平滑处理
• 边缘点的检测
• 边缘定位

边缘检测的方法：

• 基于梯度算子，如prewitter，sobel
• loG算子检测（高斯拉普拉斯）又称为Marr-Hildreth检测器
• Canny算子检测

梯度算子

梯度算子用于计算梯度偏导数的滤波器模板。（又叫差分算子、边缘算子、边缘检测算子）

常用梯度算子：

• Roberts：

• Prewitt：

• Sobel:

缺点：为对图像噪声和边缘本身特性采取预防措施

LOG算子

算子为$\nabla^2 G$，为拉普拉斯算子作用在高斯算子上。形如墨西哥草帽，数字化之后的核为：

检测步骤：

• 用一个$n\times n$的高斯低通滤波器对输入图像进行滤波
• 计算有第一步得到的图像的拉普拉斯
• 找到拉普拉斯结果中的零交叉点（正负交界处，检测某点与其邻域是否有正负变化）

寻找高斯滤波器核尺寸的方法：核尺寸应为大于六倍标准差的最小奇整数（即$6\sigma$）

Canny算子

特点

• 低错误率，所有的边缘点都被找到
• 边缘点被很好的定位
• 单一的边缘响应（正式边缘周围只返回一个点）

过程：

• 用高斯滤波器平滑输入图像
• 计算梯度复读图像和角度图像
• 对复读图像应用非最大抑制
• 用双阈值处理和连接分析来检测并连接边缘

连接

边界连接的方法：

• 局部处理
• 区域处理
• 全局处理

局部处理

过程

• 计算输入图像f(x,y)的梯度幅度矩阵M(x,y)和梯度角度矩阵a(x,y)
• 根据公式形成一幅二值图像，下式A是一个指定的角度，TA表示可以接受的范围

• 扫描二值图像的行，在不超过指定长度K的每一行中填充所有的缝隙（缝隙一定要在一个1或多个1的两端，对每一行独立处理）
• 为在任何其他方向上检测缝隙，以该角度旋转该图像，应用上一步的水平扫描过程，将结果以反方向旋转回来

全局处理

Houge变换：可以用于将边缘像素连接起来得到边界曲线，优点在于受噪声和曲线间断的影响较小

在已知曲线形状的条件下，Hough变换实际上是利用分散的边缘点进行曲线逼近，它也可看成是一种聚类分析技术

用简单的例子解释：为确定空间中一条直线$y = ax+b$只需要确定两个参数a和b，但为通过数个离散的点来拟合直线，可以做如下变换：将a和b作为横纵坐标，y和x作为待定参数。

通过在新坐标中的交点来确定对应的参数a和b，继而拟合出对应的直线。而真正的霍夫变换是通过极坐标的直线表示形式进行的：即$xcos\theta +ysin\theta = \rho$

对于图像的每一个亮的像素。做过该点的直线，进行霍夫变换，得到霍夫变换图选一个阈值，图中大于该阈值的点即为找到的直线

opencv中的hough变换直线检测函数为：

void cv::HoughLines	(	InputArray	image,
                        OutputArray	lines, # 储存了霍夫线变换检测到线条的输出矢量，以（rho，theta）表示，坐标原点为左上角
                        double	rho, # 以像素为单位的距离精度。另一种形容方式是直线搜索时的进步尺寸的单位半径
                        double	theta, # 以弧度为单位的角度精度。另一种形容方式是直线搜索时的进步尺寸的单位角度。
                        int	threshold, # 累加平面的阈值参数，即识别某部分为图中的一条直线时它在累加平面中必须达到的值
                        double	srn = 0, # 进步尺寸的除数距离
                        double	stn = 0, # 进步角度的除数距离
                        double	min_theta = 0,
                        double	max_theta = CV_PI )

lines = cv2.HoughLines(	image, rho, theta, threshold, lines, srn, stn, min_theta, max_theta)

分割

分割的依据

• 不连续性：以灰度的突变为基础分割一幅图像
• 相似性：以灰度相似性将图像分为相似的区域，如区域生长、区域聚合等

阈值分割公式：

不同的阈值处理方法:

• 全局阈值处理：当T是适用于整图的常数
• 可变阈值处理：当T在一幅图像上可以改变时
• 局部阈值处理：任何点出的阈值取决于其邻域的特性（如邻域平均灰度）
• 动态阈值处理：如果T取决于空间坐标（x，y）本身，也成为自适应阈值处

影响分割的因素：

• 波峰间的间隔
• 噪音
• 物体与背景的相对尺寸
• 光源均匀性
• 图像反射特性的均匀性

全局阈值处理

基本全局阈值分割

当物体和背景灰度十分明显时可以采用，找到合适的阈值的步骤：

• 估计一个初始阈值
• 用该初始阈值分割图像，产生$G_1、G_2$
• 对$G_1、G_2$计算平均灰度值$m_1、m_2$
• 计算新阈值$T=\frac{m_1+m_2}{2}$
• 重复上述2~4步，直到迭代的差小于一个给定参数$\Delta T$

Otsu最佳全局阈值分割

步骤：

• 计算输入图像的归一化直方图，使用$p_i,i=0,1,…,L-1$表示直方图各个分量
• 计算累计和$P_1(k) = \underset{i=0}{\overset{k}{\Sigma}}p_i$$
• 计算累计均值$m(k) = \underset{i=0}{\overset{k}{\Sigma}}ip_i$
• 计算全局灰度均值$m_G = \underset{i=0}{\overset{l-1}{\Sigma}}ip_i$$
• 计算类间方差$\sigma_B^2(K) = \frac{[m_GP_1(k)-m(k)]^2}{P_1(k)[1-P_1(k)]}$
• 得到Otsu阈值k*，即使类间方差最大的k值
• 计算k*处的可分性度量$\eta(k) = \frac{\sigma_B^2(k)}{\sigma_G^2}$

区域阈值分割

区域生长

根据事先定义的准则将像素或者子区域聚合层更大区域的过程

从一组生长点开始（可以是单个像素也可以是某个小区域），将与该生长点性质相似的相邻像素或者区域与生长点合并，形成新的生长点，重复此步骤知道不能生长为止

步骤：

• 选择合适的种子点
• 确定相似性准则
• 确定生长停止条件

区域分裂与聚合

区域分裂原理

• 选择一个分裂条件Q，如果区域不满足Q，则将区域分裂层四个象限区域
• 从整图开始持续分裂，直到所有区域都满足Q为止
• 通常需要先定一个不能再进一步分裂的最小四象限区域的尺寸

聚合原理

• 对于2个相邻的区域，仅当两个区域的并集满足条件Q时才合并
• （小区域满足Q，不见得2个满足Q的小区域合并之后同样满足Q）
• 当无法进一步聚合时，停止操作

• 新章
  数据结构笔记
• 旧篇
  现代机器人学（运动学）